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NOMENCLATURE 


dip.  e) 


coefficients 

velocity;  phase  velocity 

dimensionless  imaginary  number  in  wave  potential  for  refracted  wave  when 
angle  of  refraction  is  complex;  function;  frequency 

V=T 

,  2n 
wave  number  =  — 

A 

source  term  in  full  wave  theory 
index  of  refraction;  mode  number 

ray  parameter  =  sin0/c;  time  level  in  finite  difference  scheme 

&& TfX  — 

range,  =  \/x2+y2  or  =  \fx2  +y2  +z2 ;  radial  distance  from  the  center  of  the  earth 
radius  of  the  earth 

imaginary  parameter  analogous  to  wave  frequency 


w, ie  displacement 

notation  for  displacement  at  a  point  (m  Ax,  n  Ae)  at  a  time  pAt  in  finite  differ¬ 
ence  scheme 

x,  v,  z  Cartesian  spatial  coordinates 

Ax,  Az,  At  spacing  between  points  in  a  space-time  grid  in  finite  difference  scheme 
A,  B  amplitudes 

C  velocity;  amplitude 

Gfp)  product  of  all  reflection  and  transmission  coefficients  for  a  ray 

velocity  of  Rayleigh  wave  at  the  surface  of  a  half-space 
[G]  damping  matrix  in  finite  element  scheme 

F.  amount  of  energy  per  unit  time  per  unit  solid  angle 
F(w)  source  time  function 

G  Green’s  function 

H  Heaviside  function;  width  of  a  waveguide  or  thickness  of  a  layer 

H0  Hankel  function 

/  intensity  of  wave  (amount  of  energy  per  second  crossing  a  unit  area) 

A'0  modified  Bessel  function 

[^1 1  Ml.  (G]  stiffness  matrix  or  related  matrix  in  finite-element  scheme 
M  amplitude 

[M]  inertia  matrix  in  finite-element  scheme 

P  pressure  field 

I P !  time-independent,  complex  force  amplitude 

|  Q\  column  vector  representing  external  forces  applied  to  nodal  points  of  a  finite 

element  grid 

s/r2  +zi ;  reflection  coefficient 

pp  reflectivity  (composite  reflection  coefficient) 

S  coordinate  distance  measured  along  a  ray;  amplitude 

T  period  of  wave;  phase  function;  transmission  losses;  travel  time 

U  group  velocity;  displacement 

V  velocity 

W  wave  front 

X  Z  Cartesian  spatial  coordinates 
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IV 


compressions!  wave  velocity 
shear  wave  velocity 
Lame  constant;  wavelength 
Lamd  constant  (rigidity) 

X+2/j 

density 

2tt 

angular  frequency  =—  =  Vk 

angular  separation  between  source  and  receiver  measured  from  center  of  earth 

travel  time  of  a  ray;  phase  shift 

T(p)-pX(p) 

angle  of  incidence  or  angle  of  emergence 
prA+rip) 

imaginary  part  of  complex  angle  of  refraction;  displacement 

column  vector  representing  displacement  field  of  nodal  points  in  finite-element  grid 

parameter  analogous  to  radial  wave  number;  angle  of  incidence 

^,'4™ sin ? 

time  series  in  full  wave  theory 

((l/aJ)-pJ]Vi 

arbitrary  radius 

wave  potential 

wave  potential 

[X(jc)+2p(x)]  (3m/3x),  P- waves  or  p(x)  (du/dx),  5-waves 

3.i4159... 

convolution 

imaginary  part  of  a  complex  number  or  function 
real  part  of  a  complex  number  or  function 


REVIEW  OF  METHODS  FOR  GENERATING 
SYNTHETIC  SEISMOGRAMS 


Lindamae  Peck 


SECTION  1.  INTRODUCTION 

A  seismogram  is  a  time  record  of  ground  motion.  Displacements  from  the  rest  position  of  the 
recording  instrument  correlate  with  the  arrival  of  various  seismic  waves  at  the  location  of  the  in¬ 
strument,  although  identification  of  the  waves  and  determination  of  precise  arrival  times  may  be 
difficult  because  of  the  presence  of  noise  on  the  seismogram.  The  characteristics  of  a  seismogram 
(wave  arrivals,  time  between  arrivals,  amplitude  and  direction  of  displacements)  are  determined 
by  the  nature  of  the  source,  the  structure  of  the  earth  along  the  wave  path,  and  the  characteristics 
of  the  receiver  (recording  instrument).  The  equivalent  information  is  required  to  create  a  synthet¬ 
ic  seismogram:  a  mathematical  representation  of  the  source,  a  model  of  the  geology  between 
source  and  receiver,  and  a  surface  position  designated  as  the  location  of  the  receiver,  which  is 
modeled  mathematically  in  terms  of  its  response  characteristics.  The  model  of  subsurface  geology 
is  usually  simplified  to  represent  gross  structural  features  and  ignore  details  of  the  geology  (if 
known)  in  the  vicinity  of  the  source  and  receiver.  The  characteristics  of  a  synthetic  seismogram 
are  most  sensitive  to  features  (velocity  and  density  profiles,  layers,  irregular  structure)  of  the  geo¬ 
logical  model  employed. 

The  usefulness  of  a  synthetic  seismogram  lies  in  how  closely  it  duplicates  actual  seismograms. 
The  degree  of  correlation  between  the  two  is  a  measure  of  the  accuracy  of  the  geological  model 
used  in  computing  the  synthetic  seismogram.  When  a  structure  on  the  scale  of  the  Earth  itself  is 
modeled,  an  Earth  model  that  results  in  a  close  match  between  synthetic  and  actual  seismograms 
is  unlikely  to  be  unique.  Usually,  however,  a  selection  among  the  models  can  be  made  based  upon 
independent  criteria  such  as  consistency  with  petrological  or  thermodynamic  models  of  the  Earth’s 
interior.  When  structure  on  the  scale  of  a  few  kilometers  or  less  is  modeled,  the  problem  of  the 
uniqueness  of  the  model  is  less  and  the  model  is  often  sufficiently  well  constrained  on  the  basis 
of  surface  geology  or  seismic  surveys. 

Synthetic  seismograms  are  constructed  for  several  purposes.  They  are  useful  in  investigations 
of  the  geology  of  an  area.  The  correlation  between  recorded  and  synthetic  seismograms  aids  in 
the  construction  and  refinement  of  a  geologic  model  of  the  area.  Based  on  a  preliminary  model  so 
obtained,  the  deployment  of  seismic  surveys  to  resolve  the  subsurface  structure  can  be  done  ef¬ 
ficiently  and  relatively  economically.  In  this  way,  a  synthetic  seismogram  is  useful  as  a  recon¬ 
naissance  tool  and  provides  information  necessary  for  engineering  site  studies.  Other  applications 
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of  synthetic  seismograms  require  a  previous  knowledge  of  the  geology.  When  a  sufficiently  accu¬ 
rate  geologic  model  is  known,  synthetic  seismograms  can  be  used  to  predict  the  effects  of  explo¬ 
sions  or  earthquakes  (the  amplitude  of  the  ground  motion  at  a  site  due  to  a  source  of  given  mag¬ 
nitude)  or  to  derive  a  theoretical  model  of  vehicle-generated  or  artillery-generated  ground  motion 
for  use  in  detection  and  classifier  design. 

The  validity  of  a  synthetic  seismogram  is  dependent  upon  the  accuracy  of  the  geologic  model 
employed,  the  applicability  of  the  assumptions  inherent  in  the  derivation  of  the  wave  equation, 
and  the  approximations  incorporated  in  the  solution  to  the  wave  equation. 

The  wave  equation  is  derived  from  a  statement  of  the  balance  of  forces  acting  on  a  unit  volume 
of  unbounded  medium  (e.g.  Grant  and  West  1965,  Officer  1974,  Aki  and  Richards  1980).  Several 
assumptions  generally  are  made  that  result  in  a  wave  equation  of  the  form 


a?=c''(a?+a72+a?) 


where  r  is  time,  c  is  velocity,  and  jr,r  and  z  are  spatial  (Cartesian)  coordinates.  One  assumption  in 
deriving  this  equation  is  that  external  forces  such  as  gravity  acting  on  the  medium  may  be  neglect¬ 
ed.  A  second  is  that  the  medium  is  linearly  elastic;  this  determines  the  form  of  the  stress-strain 
relation  that  is  substituted  in  the  balance  of  force  equation.  A  third  is  that  the  medium  is  homo¬ 
geneous  and  isotropic  and  therefore  the  elastic  parameters  (e.g.  the  Lame'  constants)  of  the  medium 
are  not  functions  of  the  spatial  coordinates. 

A  synthetic  seismogram  is  generated  by  solving  the  wave  equation  subject  to  the  constraint 
that  the  boundary  conditions  of  the  problem  are  satisfied.  The  boundary  conditions  pertain  to 
the  components  of  stress  and  displacement  at  the  free  surface  of  the  medium,  at  any  interfaces 
within  the  medium,  and  at  infinite  depth  in  the  medium.  Various  analytic  and  numerical  methods 
are  available  for  use  in  solving  the  wave  equation;  a  distinction  among  them  is  that,  in  general, 
analytic  methods  yield  an  approximate  solution  to  the  wave  equation  while  numerical  methods 
yield  an  exact  solution  to  an  approximation  of  the  wave  equation.  The  selection  of  a  particular 
method  is  based  upon  the  wave  types  that  are  to  be  represented  in  the  seismogram  and  on  the  de¬ 
gree  of  structural  complexity  that  the  method  of  generating  the  seismograms  is  capable  of  han¬ 
dling. 

This  report  presents  a  description  and  evaluation  of  commonly  used  methods  of  generating 
synthetic  seismograms.  Us  intent  is  to  indicate  the  strengths  and  weaknesses  of  the  methods  and 
to  provide  examples  of  recent  applications  of  them.  The  cited  literature  serves  as  an  initial  refer¬ 
ence  for  individuals  to  use  in  pursuing  an  approach  to  a  particular  problem.  A  concern  of  this  re¬ 
port  is  to  identify  methods  capable  of  generating  synthetic  seismograms  for  physical  situations 
such  as  thin  beds  or  near-surface  velocity  inversions  that  are  typical  of  winter  seismological  condi¬ 
tions.  The  mathematical  representation  of  source  characteristics  and  receiver  response  is  not  con¬ 
sidered  here;  the  cited  references  and  general  texts  such  as  that  of  Aki  and  Richards  (1980)  show 
how  such  factors  are  incorporated  in  the  process  of  generating  synthetic  seismograms. 


SECTION  2.  WAVE  PROPAGATION  IN  THE  EARTH 

The  types  of  seismic  waves  are  distinguished  by  particle  motion  and  velocity  of  propagation. 
The  compressional  wave  (or  synonymously  the  irrotational  wave,  longitudinal  wave,  dilatational 
wave  or  P- wave)  propagates  at  a  velocity  of  a  =  \/(X+2p)/p  through  an  isotropic,  homogeneous 
medium  that  is  alternately  condensed  and  rarefied;  the  particle  motion  is  parallel  to  the  direction 
of  wave  propagation.  Here,  X  and  p  are  the  Lame  constants  and  p  is  the  density  of  the  medium. 
The  shear  wave  (or  synonymously  the  rotational  wave,  transverse  wave  or  S-wave)  propagates  at 


a 


a  velocity  (3  =  y/pip:  distortion  of  the  medium  is  in  the  plane  perpendicular  to  the  direction  of 
propagation  and  is  designated  as  SH  (horizontal)  or  S  V  (vertical)  depending  on  its  orientation. 

The  P-  and  5-waves  are  body  waves. 

Rayleigh  waves  are  coupled  P-SV  waves  with  motion  that  at  the  surface  defines  an  ellipse  per¬ 
pendicular  to  the  surface  and  that  exponentially  decays  in  amplitude  with  depth;  the  waves  are 
dispersive  (propagation  velocity  depends  on  frequency)  in  a  layered  medium  and  nondispersive 
(with  a  velocity  Cr  equal  to  0.9194/3  when  Poisson’s  ratio  is  14)  when  propagating  at  the  surface 
of  a  half-space  of  shear  velocity  /?.  Love  waves  are  57/  waves  that  propagate  sinusoidally  within  a 
layer,  decay  exponentially  with  depth  outside  the  layer,  and  exhibit  dispersion;  the  velocity  of 
the  waves  approaches  the  shear  velocity  of  the  surface  layer  in  the  limit  of  high  frequencies  and 
approaches  the  (higher)  shear  velocity  of  the  underlying  half-space  in  the  limit  of  low  frequencies. 
Rayleigh  waves  and  Love  waves  are  called  surface  waves.  They  are  generated  along  with  body 
waves  (unless  the  source  is  located  at  a  depth  such  that  surface  waves  are  not  excited)  and  propa¬ 
gate  horizontally  at  the  surface  with  an  attenuation  due  to  geometrical  spreading  that  is  propor¬ 
tional  to  l/(range)'^.  Since  body  waves  are  attenuated  as  \/r  and  head  waves  (see  below)  as  1/r 2 , 
the  major  portion  of  the  wave  energy  reaching  points  distant  from  the  source  is  carried  in  surface 
waves.  Surface  waves  are  viewed  as  noise  to  be  eliminated  when  vertically  or  near  vertically  inci¬ 
dent  (reflected)  waves  are  of  interest,  as  is  typically  the  case  in  exploration  seismology,  or  as  use¬ 
ful  counterparts  to  body  waves  when  analyzed  to  extract  information  on  source  parameters  and 
earth  structure  between  source  and  receiver.  Because  of  the  conditions  under  which  they  exist 
and  the  manner  in  which  they  propagate,  surface  waves  are  exceptionally  well  suited  for  use  in 
the  delineation  of  near-surface  layered  structure. 

Diffracted  waves  are  present  when  plane  waves  are  incident  on  a  feature  that  has  a  radius  of 
curvature  less  than  or  equal  to  the  wavelength  of  the  incident  waves.  A  discontinuous  boundary 
(such  as  a  geologic  bed  displaced  by  a  fault)  or  a  nonplanar  boundary  commonly  causes  diffracted 
waves  to  be  present  in  the  wavefield.  At  distances  from  the  diffracting  source  that  are  large  rela¬ 
tive  to  wavelength,  Huygens'  construction  (e.g.  Halliday  and  Resnick  1974)  is  used  to  determine 
the  diffracted  wave  energy  at  a  particular  location. 

A  point  source  in  a  homogeneous,  isotropic  medium  generates  /’-waves  (and  5-waves  if  the 
medium  is  nonfluid)  that  radiate  in  spherical  fronts;  i.e.  the  waves  propagate  in  all  directions 
from  the  source  and  points  of  the  medium  that  are  in  the  same  phase  of  motion  are  located  on 
surfaces  (wave  fronts)  that  define  spheres.  As  a  wave  propagates  away  from  the  source,  the  radius 
of  curvature  of  its  wave  fronts  increases  and  the  wave  fronts  become  increasingly  planar  over  a 
limited  region  (e.g.  Halliday  and  Resnick  1974,  Fig.  1).  The  wave  fronts  of  a  plane  wave  are  paral¬ 
lel  planes  perpendicular  to  the  direction  of  propagation  of  the  wave;  the  rays  are  parallel  lines  nor¬ 
mal  to  the  wave  fronts.  Figure  2  shows  wave  fronts  and  rays  of  a  plane  wave  propagating  in  the 
x-z  plane  at  a  velocity  V\  the  propagation  direction  of  the  plane  wave  is  at  an  angle  Q  to  the  z-axis. 


\ 

Figure  I.  Wave  fronts  shown  in  two 
'  dimensions  as  dashed  lines  emanating 
I  from  the  point  source,  S.  Radial 
arrows  normal  to  the  wave  fronts  are 
/  rays.  Each  ray  indicates  the  direction 
of  propagation  of  the  wave.  The  spheri¬ 
cal  wave  fronts  become  locally  planar 
as  distance  from  S  increases. 
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Figure  2.  Wave  fronts  of  a  plane  wave  shown 
as  dashed  lines.  Rays  are  shown  as  arrows. 


The  plane  wave  may  be  represented  by  the  potential 

O  =  .1  exp  Hkx. v  +  kyz  -  re/)  (1) 

where  .1  is  the  amplitude  of  the  wave  potential:  Ax  =  k  sinO  and  k  =  k  cos 0  are  the  horizontal  and 
vertical  components  of  the  wave  number  (where  k  =  2rr/A  and  X  is  wavelength),  w  =  2ttlT  =  Vk  is 
angular  frequency  and  T  is  period. 

In  a  homogeneous,  isotropic,  unbounded  medium,  the  only  waves  that  arrive  at  the  receiver 
iiom  an  isotropic  source  are  direct  waves  (Tig.  3a).  In  a  homogeneous,  isotropic  half-space,  surface- 
reflected  waves  (and  Rayleigh  w  aves)  also  arrive  at  the  receiver  (Fig.  3b).  If  velocities  (a,  0)  in¬ 
crease  uniformly  w  ith  depth  in  the  half-space,  refracted  waves  arrive  at  the  receiver  after  traveling 
minimum-time  paths  (Fig.  3c).  If  the  medium  consists  of  homogeneous  layers  of  contrasting  den¬ 
sity  and/or  velocity ,  waves  incident  at  an  interface  generate  reflected  waves  and  transmitted  waves; 
which  waves  arrive  at  the  receiver  is  determined  by  the  angle  of  incidence  and  the  number  of  times 
the  waves  are  reflected  (Fig.  3d). 

A  /'-wave  or  a  .S' I '-wave  that  is  incident  obliquely  at  a  planar  boundary  between  two  media  in 
welded  contact  generates  two  reflected  waves  (/',  ST)  and  two  transmitted  waves  (P,  SI).  The 
angle  at  which  each  reflected  or  transmitted  wave  leaves  the  interface  is  determined  by  a  general¬ 
ized  form  of  Snell’s  law.  Figure  4  shows  a  P- wave  ray  in  the  x-z  plane  arriving  at  a  planar  bound¬ 
ary  with  an  angle  of  incidence,  0F  .  By  Snell’s  law, 
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where  inc,  rfl  and  rfr  denote  incident,  reflected  and  refracted  (transmitted)  waves,  respectively. 
The  constant  is  termed  the  ray  parameter  p.  If  a  SC-wave  in  the  x-z  plane  is  incident  at  an  angle 
0S  ,  the  governing  equation  is 
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A  S//- wave  in  the  v-r  plane  incident  at  the  boundary  does  not  generate  P- waves  or  SC-waves  be¬ 
cause  its  motion  is  parallel  to  the  boundary.  Equations  similar  to  cq  1  may  be  written  for  the 
cornpressional  wave  potential  and  the  shear  wave  potential  in  each  medium.  The  potentials  are 
substituted  in  equations  stating  the  continuity  of  stress  and  of  displacement  across  the  boundary; 
the  set  of  equations  is  then  solved  for  the  coefficients  of  reflection  and  of  transmission  (refrac¬ 
tion).  Each  coef  ficient  is  a  ratio  of  amplitudes  (of  two  potential  equations)  and  depends  upon 
the  angle  of  incidence  of  the  waves  and  on  the  elastic  properties  and  densities  of  the  media.  Co¬ 
efficients  may  be  obtained  also  in  terms  of  the  amplitudes  of  displacement  or  energy.  Graphs  of 
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Figure  3.  Wave  types  ci  bounded  and  un¬ 
bounded  media.  S  denotes  source,  K  de¬ 
notes  receiver. 

a.  Direct  wave  in  infinite, 
homogeneous  medium. 

b.  Surface-reflected  wave, 
e.  Refracted  wares. 

d.  Multiply  reflected  waves 
in  a  layered  medium. 


Figure  4.  Ray  for  each  reflected  /P,  S )  plane 
wave  and  each  transmitted  ( P,  S)  plane  wave 
generated  by  a  plane  P-trare  incident  at  an 
angle  0Pint  at  a  planar  boundary.  The  P- 
and  S -wave  velocities  are  a,  and  in  med¬ 
ium  1  and  ot2  and  $2  in  medium  2. 


transmission  and  reflection  coefficients  as  a  function  of  angle  of  incidence  (for  various  density  and 
velocity  contrasts)  are  given  in  Grant  and  West  (1965,  pp.  57-59)  and  Telford  et  al.  (1978,  p.  254). 
As  the  angle  of  incidence  of  the  P- wave  increases  from  0 .>  =  0,  the  angle  of  refraction  of  the 

'inc  °  0 

transmitted /’•wave  also  increases  (eq  2).  When  0,,  =  sin  (a,/a2),  0P  ,  =90  and  the  trans- 

mitted  P- wave  propagates  in  medium  2  at  grazing  incidence;  0F.^  is  termed  the  critical  angle 
(0„  )  for  refraction  of  the  P-wave  and  the  refracted  wave  is  termed  the  head  wave.  The  analysis 

'  cr it 

of  head  waves  is  based  upon  the  theory  ot  curved  wave  fronts  at  a  planar  interface  (or  plane  waves 
incident  at  curved  interfaces).  The  reflected  /’-wave  energy  is  maximum  at  9,,  =9.,  .  because 

this  is  the  condition  for  total  internal  reflection. 

When  (i|,inc>0t,crlt-  Snell's  law  would  require  (a2/a1)sinOP.nc  be  greater  than  1.  In  this  case, 
0(>  fr  is  taken  to  be  complex,  0\>rfr  =  n/2  -  ib  and  b  =  cosh'1  [(a2/a,)sinf)pjnc]  (Grant  and  West 
ll>65.  p.  60).  When  0|'rfr  is  substituted  in  a  potential  equation  for  the  transmitted  P-wave,  the 
equation  that  results  has  the  form 

=  //exp  [/  W  sin 0,,  (x+r/)|  . 
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w  here 
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is  an  imaginary  number.  Therefore,  the  transmitted  /'-wave  as  expressed  by  is  not  an  ordi¬ 
nal  /'  wave,  but  rather  is  an  inhomogeneous  wave  characterized  by  an  exponential  decrease  in 
amplitude  with  vertical  distance  from  the  interface.  In  terms  of  the  ray  parameter,  the  transmit¬ 
ted  /’-wave  is  an  inhomogeneous  wave  when  p  >  1  fa2  .  If  02  <  Qi ,  a  P- wave  incident  at  an  angle 
of  dp  >  /),,  generates  a  transmitted  .ST-wave  at  a  real  angle  of  refraction  and  so  the  .5 1’-wave 

T  inc  '  cr  it  - 

propagates  away  from  the  boundary  .  If  02  >at ,  there  is  a  critical  angle  for  refraction  of  the  .ST- 
wave,  When  Of  =  0SV  .  =  sin  1  (a,  P2  ).  the  .ST  -wave  generated  by  the  incident  /’-wave  is  a 

he  id  wive  Ilf)','."'  "’f)sv  .  the  angle  of  refraction  of  the  .ST-wave  is  complex  and  the  trans- 
mitted  .ST-wave  is  an  inhomogeneous  wave.  Both  /’and  .ST  transmitted  motion  are  trapped  at  the 
interlace. 
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Figure  11.  Plane  wave  representation  of  normal  modes  of 
the  first  three  orders.  The  boundary  conditions  are  a  pres¬ 
sure-released  surface  and  a  pressure-reinforcement  (rigid) 
bottom  (adapted  from  Urick  1979). 

decay  with  range.  If  the  imaginary  part  of  the  complex  wave  number  is  much  smaller  than  the 
real  part,  then  the  real  part  is  used  in  determining  phase  velocity  (C  =  T/sin7n  =  co/An)  and  group 
velocity  (U  =  do:/dkn ),  while  the  imaginary  part  is  a  damping  coefficient  (Brekhovskikh  1960,  p. 
405 ).  Reflection  coefficients  that  are  dependent  on  y  are  used  when  the  waveguide  boundaries  are 
not  perfect  reflectors. 

The  approach  used  above  to  determine  the  condition  for  constructive  interference  of  plane 
waves  in  a  homogeneous  waveguide  can  be  applied  in  the  case  of  a  layered  waveguide  (velocity  de¬ 
pendent  on  depth  within  waveguide).  This  requires  that  the  effect  of  partial  reflection  at  the  lay¬ 
er  interfaces  on  the  phase  of  each  plane  wave  and  on  the  energy  content  of  the  plane  waves  be  in¬ 
corporated  in  the  solution. 

Normal  modes  that  have  S/1  polarization  are  Love  waves,  which  are  specified  by  their  mode 
number  (ti  =  0,  1.2,  ...).  The  fundamental  normal  mode  with  N'T  polarization  corresponds  to  the 
dispersive  Rayleigh  wave  and  is  designed  as  the  Rayleigh  mode;  the  higher  order  modes  are  desig¬ 
nated  as  the  first,  second,  etc.,  NT  mode  (Grant  and  West  1965). 

There  are  many  computer  programs  available  to  obtain  normal  mode  solutions  to  underwater 
wave  propagation  problems,  some  of  which  are  useful  in  solid  earth  studies.  The  normal  mode 
solution  can  he  obtained  by  solving  the  Helmholtz,  equation;  this  is  a  time-independent  wave  equa¬ 
tion  that  results  when  the  pressure  field  is  assumed  to  be  of  constant  angular  frequency,  P  =  S(r,z) 
exp  (-ieot).  If  the  spatial  factor  has  the  form  S  =  //(/'  (kQr)  f(r,  r)  and  the  asymptotic  approxima¬ 
tion  of  the  Hankel  function  (valid  for  A0r»  I )  is  substituted,  then  the  Helmholtz  equation  that 
results  (in  cylindrical  coordinates)  is 
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(nodes),  the  lines  of  positive  signs  are  wave  fronts  of  positive  wave  amplitude,  and  the  lines  of 
negative  signs  are  wave  fronts  of  negative  amplitude.  Each  wave  front  impinges  on  the  waveguide 
boundaries  at  an  angle  y.  In  order  for  superimposed  plane  waves  to  interfere  constructively,  the 
phase  difference  between  the  waves  must  be  an  integral  multiple  of  2n  (2nn,n  =  0,  1 , 2, This 
requirement  leads  to  a  relationship  between  thickness  of  the  waveguide,  angle  of  incidence  of  tire 
plane  waves,  and  frequency.  Only  waves  of  certain  frequencies  will  interfere  constructively  and 
give  rise  to  a  standing  wave  pattern  in  the  z -direction;  frequency  is  uniquely  determined  when 
mode  number  (/;)  and  y  are  specified.  Tire  apparent  velocity  (velocity  in  the  .v -direction)  of  the 
plane  waves  is  C  =  F/sinT.  this  relationship  is  used  to  rewrite  the  equation  for  constructive  inter¬ 
ference  in  terms  of  C  instead  of  y.  The  result  is  a  dispersion  equation  (phase  velocity-frequency 
relation)  for  each  normal  mode. 

Since  a  critical  angle  of  incidence  exists  for  the  retention  of  wave  energy  in  the  waveguide  by 
total  internal  reflection,  there  is  a  cutoff  frequency  for  each  normal  mode.  At  frequencies  below 
the  cutoff  frequency,  the  normal  modes  are  attenuated  as  they  prop  ‘gate  even  in  a  lossless  medium 
due  to  the  transmission  of  wave  energy  across  the  waveguide  boundary  (leaky  modes).  The  energy 
lost  from  the  waveguide  as  transmitted  waves  may  be  returned  by  reflections  at  depth  and  so,  over 
short  ranges,  may  contribute  to  the  wave  field  at  the  receiver. 

This  treatment  of  normal  modes  requires  that  the  wave  energy  has  propagated  sufficiently  far 
from  the  source  for  the  representation  of  normal  modes  as  interfering  plane  waves  to  be  valid.  The 
reflection  coefficient  for  spherical  waves  is  the  sum  of  the  plane-wave  reflection  coefficient  and  a 
correction  term  that  goes  to  zero  at  very  high  frequencies.  The  correction  term  is  negligible  when 
source  and  receiver  are  located  at  distances  from  the  waveguide  boundary  that  are  large  relative  to 
wavelength  (Brekhovskikh  1980,  p.  247;  Clay  and  Medwin  1977,  p.  64).  In  a  mathematical  deriva¬ 
tion  of  the  normal  mode  solution,  the  asymptotic  representation  of  the  Hankel  function  is  substi¬ 
tuted  for  the  function  itself,  which  is  valid  for  distances  from  the  source  that  are  large  relative  to 
wavelength. 

The  normal  mode  solution  for  the  relatively  simple  case  of  a  fluid  layer  over  a  rigid  half-space 
(Fig.  10)  is  given  here  as  an  example  of  the  solution  format.  (The  case  of  a  solid  layer  over  a  rigid 
bottom  is  more  complex  because  two  types  of  waves  are  present  \P,  5]  and  so  two  potentials  are 
required  to  express  the  normal  mode  solution.  If  the  half-space  is  nonrigid,  a  solution  must  be 
found  for  that  region  as  well,  e.g.  that  of  Ewing  et  al.  1957,  p.  138.)  The  potential  0  for  horizon- 
tails  propagating  waves  in  the  layer  is 

0  =  21  fln  [/?,  (yn)exp(-fAz  cosyn)  +  exp(/*z  cosyn)|  exp(/*.r  sinyn )  (15) 

( Brekhovskikh  1960,  p.  407),  where  n  designates  a  particular  normal  mode,  R ,  is  the  reflection  co¬ 
efficient  for  waves  incident  from  above  at  the  lower  boundary  of  the  waveguide,  and  the  coefficient 
Bn,  which  is  a  function  of  the  position  and  strength  of  the  source,  determines  the  degree  of  excita¬ 
tion  of  each  normal  mode.  The  z -dependence  of  the  amplitude  of  each  normal  mode  is  determined 
by  the  expression  in  the  brackets  of  eq  13.  The  potential  0  is  the  sum  of  potentials  lor  the  upgoing 
(0_)  and  downgoing  (0+)  plane  waves  where  0.  =  A  exp[/  (-Azz  +  k%x  )] ,  0+  =  B  cxp[/  (Azz  +  A**)] , 
Ax  =  A  siny.  and  k t  -  A  cosy.  An  equivalent  equation  could  be  written  in  terms  of  the  coefficient 
.1  and  the  reflection  coefficient  for  waves  incident  from  below  at  the  upper  boundary  of  the  wave¬ 
guide,  /?2.  The  potential  0  (eq  13)  satisfies  the  equation /?  ^2  exp  (2rAz//)  =  1 .  Lxamplesof 
plane  waves  at  three  angles  of  incidence  and  the  normal  modes  that  result  from  constructive  inter¬ 
ference  are  given  in  Figure  1  1 .  For  a  given  frequency,  a  higher  order  mode  corresponds  to  a  steep¬ 
er  ray  path. 

The  normal  mode  solution  to  the  wave  equation  can  be  generalized  to  include  attenuation  due 
to  propagation  in  an  imperfectly  elastic  (lossy )  medium  and  imperfect  reflections  at  the  waveguide 
boundaries;  this  is  done  by  making  the  wave  number  complex  and  so  introducing  an  exponential 


Consider  a  homogeneous  elastic  layer  of  thickness  fl  over  a  homogeneous  half-space,  each  with 
shear  wave  velocity  (/?)  and  density  (p)  as  shown  in  Figure  8.  The  dispersion  equation  for  Love 
waves  is 
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its  solution  is  graphed  in  Figure  9.  The  cutoff  frequency  is  zero  for  the  fundamental  mode  (/;  =  0) 
and 
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tor  the  higlier  modes  (n  >  0)  ( Aki  and  Richards  1980,  p.  264).  For  a  given  mode,  the  phase  veloc¬ 
ity  is  equal  to  (?2  at  the  cutoff  frequency  and  is  equal  to  &  in  the  limit  of  or*00  (asymptotic  solu¬ 
tion  to  eq  12 ),  This  is  physically  reasonable  in  that  high  frequency  waves  (short  wavelengths)  are 
confined  to  the  vicinity  of  the  layer  surface  whereas  low  frequency  waves  (long  wavelengths)  pene¬ 
trate  the  deeper,  high  velocity  medium  and  are  little  affected  by  a  relatively  thin  surface  layer.  A 
dispersion  relation  for  Rayleigh  waves  is  more  complicated.  Pilant  (1979),  Takeuchi  and  Saito 
(1972),  and  Aki  and  Richards  (1980)  treat  Rayleigh  wave  propagation  in  a  layered  medium.  Equa¬ 
tions  11  a,  b  indicate  that  Rayleigh  waves  are  retrograde  elliptical  at  the  surface  and  prograde  ellipti 
cal  below  the  nodal  plane.  Results  of  computer  calculations  of  surface  wave  dispersion  indicate 
that  higher  order  modes  are  prograde  elliptical  at  the  surface. 

The  characterization  of  Rayleigh  waves  and  Love  waves  in  terms  of  modes  indicates  a  relation¬ 
ship  between  surface  waves  and  normal  modes;  the  relationship  will  become  apparent  in  the  discus¬ 
sion  of  normal  mode  theory  in  the  next  section.  Here,  it  should  be  noted  that  certain  conditions 
must  be  satisfied  if  the  approach  of  modeling  surface  waves  as  plane  waves  that  are  multiply  re¬ 
flected  in  a  layer  is  to  be  valid.  The  use  of  plane  waves  requires  that  the  waves  have  propagated 
sufficiently  far  from  the  source  to  be  well  represented  as  planar  rather  than  spherical  or  cylindrical. 
This  condition  is  expressed  in  terms  of  the  layer  thickness:  the  distance  from  the  source  must  be 
large  in  terms  of  the  layer  thickness  (Grant  and  West  1965,  p.  77). 


SECTION  5.  NORMAL  MODES 


A  discussion  of  normal  modes  is  given  in  Ewing  et  al.  (1957),  Grant  and  West  (1965),  Officer 
( 1974),  and  Brekhovskikh  (1980).  The  propagation  of  normal  modes  can  be  represented  as  arising 
from  the  interference  of  plane  waves  that  are  multiply  reflected  at  the  boundaries  of  the  waveguide 
One  set  of  reflected  waves  in  an  elastic  waveguide  with  plane,  parallel  boundaries  overlying  a  rigid 
half-space  is  shown  in  Figure  10.  The  solid  lines  are  wave  fronts  of  zero  amplitude  of  the  waves 
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f  igure  III.  Wave  [runts  impinging  on  boundaries  of  waveguide  at  angle  y 
(adapted  from  Crick  IVS.i). 
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(u)  group  velocity 
( c)  phase  velocity 


Figure  7.  Dispersion  curves  for  three  modes  ( I,  2,  3).  The 
cutoff  frequencies  are  cjc1,  wCJ,  and  respectively. 

Die  curves  are  drawn  for  the  situation  of  increasing  seismic 
velocity  with  depth  (adapted  from  Officer  1974). 
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Figure  9.  Graphic  solution  of  eq  12  for  the 
dispersion  of  Love  waves  in  a  single  layer 
over  a  half-space.  Each  side  of  eq  12  is  plot¬ 
ted  as  a  function  ofH[l-(Pi2/C2)]'/z/(li 
(adapted  from  Akiand  Richards  1980). 


Figure  8.  Geologic  model  for  which  the  dis¬ 
persion  equation  for  Love  waives  (eq  12)  is 
calculated. 


Surface  waves  propagating  in  a  layer  can 
be  viewed  as  arising  from  the  interference 
of  plane  waves  that  are  multiply  reflected 
at  the  boundaries  of  the  layer.  The  condi¬ 
tion  for  constructive  interference  of  the 
plane  waves  leads  to  the  development  of  a 
dispersion  equation  that  expresses  the  fre¬ 
quency-velocity  relation  of  solutions  to  the 

wave  equation  in  terms  of  the  thickness  of  the  layer,  the  velocities  and  densities  of  the  layer  and 
substratum,  and  the  angle  of  incidence  of  the  plane  waves.  The  dispersion  equation  is  multi¬ 
valued  due  to  a  tangent  function;  one  dispersion  curve  is  obtained  for  each  mode  (value  of  n, 
where  n  =  0, 1 , 2, ...,  each  corresponding  to  a  range  of  the  tangent  function)  specified.  There  is 
a  cutoff  frequency  wcnfor  all  but  the  fundamental  («  =  0)  Love  mode  and  for  all  modes  in  a  lay¬ 
ered  acoustic  medium.  At  frequencies  equal  to  or  greater  than  the  cutoff  frequency,  surface  waves 
propagate  horizontally  with  essentially  no  attenuation  in  a  lossless  medium  but  do  decay  exponen¬ 
tially  with  depth;  at  frequencies  below  the  cutoff  frequency,  the  wave  attenuates  as  it  propagates 
even  in  a  lossless  medium  due  to  the  radiation  of  energy  across  the  layer  boundary  (leaky  modes) 
as  well  as  decaying  exponentially  with  depth.  The  cutoff  frequency  increases  with  mode  number 
so  there  is  a  finite  number  of  modes  possible  at  a  given  frequency.  The  phase  velocity  (c)  is  the 
velocity  at  which  a  point  of  a  particular  wave  is  propagating.  The  group  velocity  (U)  is  the  veloc¬ 
ity  at  which  the  energy  associated  with  a  particular  frequency  group  travels.  Group  and  phase 
velocity  are  related  by  U-c-  Mdc/dX).  A  generalized  example  of  dispersion  curves  for  phase 
velocity  and  group  velocity  is  given  in  Figure  7. 
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and 


\p  =  B  exp/co  [t  -  (x/c)  +  bz\ .  (9) 

The  amplitude  coefficient  B  must  be  directed  in  the  y -direction,  so  eq  9  is  rewritten  as 

\py  =Z?exp/c.'[/  ~(x/c)  +  bz].  (10) 

Substitution  of  eq  8  and  10  in  the  wave  equation  (32 /dr2  =  u2(32/3*2 ),  v  =  a  or  (3)  leads  to  = 
1/a2  -  1/c2  and  b2  =  1//32  -  1/c2.  If  c  >  a  >  P,  a  and  b  are  real  and  0  and  0y  represent  two  waves 
propagating  in  oblique  directions  that  do  not  decay  with  depth.  In  order  to  obtain  a  representa¬ 
tion  for  waves  that  are  restricted  to  the  vicinity  of  the  surface,  it  must  be  that  c  <  a  and  c  <  /3 
so  that  a  and  b  are  imaginary.  By  choosing  the  imaginary  parts  of  a  and  b  to  be  positive  (a  = 
is/\/c2-  1/a2 ,  b  =  is/X/c^  - 1  //32 ),  the  waves  will  decay  exponentially  with  depth  as  they  propagate. 
The  system  of  equations  that  results  from  application  of  the  boundary  conditions  a22  =  oxl  =  0 
is  solved  for  c;  for  a  Poisson’s  solid  (medium  with  Poisson’s  ratio  of  !4),  c  =  C,  =  0.9194/3.  The  dis¬ 
placement  vector  u  for  Rayleigh  waves  has  scalar  components: 

ux  =4/ exp/co  (/-^(e-°-847Swz/c-  0.5773e-°-3933“2/c)  (11a) 

uz  =  -  iM  exp/co  (/  •  (-0.8475e-°-847Sw2/c  +  1 .4679e-°-3933w2/c)  .  (lib) 

These  equations  show  that  a  surface  wave  of  low  frequency  decays  more  slowly  than  one  of  high 
frequency.  At  the  surface,  ujux  =  1.468  and  at  a  depth  of  0.193X,  ux  =  0  (Kolsky  1953,  p.  20- 
21 ).  A  specific  Poisson’s  ratio  has  been  substituted  in  the  general  equations  in  order  to  make  the 
equations  given  here  more  tractable;  alternately,  Grant  and  West  (1965,  Fig.  3-10)  show  ujux\ 
and  the  depth  to  the  nodal  layer  as  a  function  of  Poisson’s  ratio.  'z  0 

In  a  similar  manner,  the  displacements  due  to  a  Love  wave  propagating  in  an  elastic  layer  over 
a  half-space  are  uy  =  A  exp/co  ( t  -  (x/c)  +  bxz)  +B  exp/co  (t  -  (x/c)  -b^z)  in  the  layer  and  uv  = 

C  exp/co  (t  -  (x/c)  +  b2z)  in  the  half-space,  with  bi  =  \/l//3i2  -  1/c2  and  b2  =  b/l/c2  -  l//322  .  The 
boundary  conditions  are  continuity  of  stress  and  displacement  at  the  layer/half-space  interface  and 
the  vanishing  of  the  stress  at  the  free  surface.  The  equation  that  results  upon  elimination  of  A,  B, 
and  C  is  the  dispersion  equation,  expressing  the  dependence  of  velocity  of  propagation  on  frequency. 

The  treatment  of  surface  waves  can  be  generalized  to  the  case  of  several  layers  over  a  half-space. 

In  each  layer  there  are  upward  and  downward  propagating  (and  decaying)  waves  that  interfere; 
only  downward  propagating  (and  decaying)  waves  exist  in  the  half  space.  Long  wavelength  waves 
penetrate  the  medium  more  deeply  and  so  have  a  phase  velocity  that  is  an  average  of  the  medium 
properties  to  a  greater  depth.  Boundary  conditions  on  motion  within  the  vertically  inhomogene¬ 
ous  medium  are  continuity  of  stress  and  displacement  across  the  layer  interfaces.  In  order  to  con¬ 
struct  a  dispersion  equation  for  the  multi-layered  medium,  the  motion  at  a  subsurface  layer  must 
be  related  to  motion  at  the  surface;  the  equations  expressing  this  relation  can  be  solved  by  numer¬ 
ical  integration  (e.g.  Runge-Kutta  method)  or  by  matrix  techniques  (e.g.  Thomson-Haskell  method) 
(Aki  and  Richards  1980,  p.  270).  Phase  velocity -frequency  curves  are  constructed  from  the  dis¬ 
persion  relation;  a  comparison  of  the  curves  and  field  data  indicates  how  accurately  the  model 
used  to  construct  the  curves  represents  the  velocity  profile  and  near-surface  structure  of  the  earth. 

Hermann  ( 1 978)  gives  annotated  computer  programs  for  surface  wave  seismology.  A  combina¬ 
tion  of  programs  supplied  will  1 )  solve  for  phase  and  group  velocities  as  a  function  of  period  for 
Love  and  Rayleigh  waves  in  a  multi-layered  medium,  2)  calculate  the  Rayleigh  wave  and  Love  wave 
displacements  in  the  layers,  and  3)  generate  synthetic  seismograms. 
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the  computational  advantages  of  evaluating  the  frequency  integral  first  and  then  integrating  over 
real  values  of  the  ray  parameter  (or  slowness).  Using  Chapman’s  approach,  Dey-Sarkar  and  Chap¬ 
man  ( 1978)  obtain  the  following  displacement  equation  valid  for  (ojprA  »  1 ): 


u(r.  A,  r)  ~ 


m(t) 


n(  2  rs  sinA),/! 


*  Im[A(r)  *  Z, 


Pv,C(p )  p,  1 

t=0  4 na\(p%q%prqt)'h  Qr  \0'\ 


where  r  =  radius 

r,  s  =  receiver  radius  and  the  source  radius,  respectively 
m(t)  =  the  source  term 

C(p)  =  the  product  of  all  reflection  and  transmission  coefficients  defining  the  ray 
A (t)  =  a  time  series 

q  =  a  function  of  ray  parameter  p  and  radius 
6(p )  =  prA+r(p) 

rip)  =  the  intercept  function  obtained  from  a  graph  of  travel  time  vs  range 

A  =  the  angular  separation  between  source  and  receiver  measured  from  the  center  of 
the  earth. 


The  summation  is  over  all  real  solutions  of  t  =  6(p). 

Choy  et  al.  (1980)  compare  synthetic  seismograms  generated  by  the  reflectivity  method  and 
by  full  wave  theory.  An  advantage  of  the  full  wave  theory  method  is  that  no  truncation  phase  is 
present  in  the  seismograms,  while  integration  over  a  limited  range  of  angles  of  incidence  causes  a 
truncation  phase  to  be  present  in  the  reflectivity  seismograms;  a  disadvantage  is  that  the  rays 
must  be  specified  in  order  to  be  included  in  the  full  wave  theory  calculations  (unlike  the  reflec¬ 
tivity  method  in  which  all  internal  multiples  and  conversions  are  automatically  included). 


SECTION  4.  SURFACE  WAVES 

Surface  waves  exist  in  the  near-surface  region  of  a  medium  having  a  free  surface.  They  decay 
exponentially  with  depth  in  the  medium  but  attenuate  laterally  with  a  1  /(range )Vl  dependence. 
The  types  of  surface  waves  are  distinguished  by  differences  in  particle  motion  and  in  conditions 
for  occurrence.  Love  waves  are  SH- waves  that  propagate  in  a  layered  medium.  Rayleigh  waves 
are  coupled  P-SV  waves  that  propagate  in  a  layered  medium  and  at  the  surface  of  a  half-space. 

In  the  case  of  a  single  layer  over  a  half-space,  Love  waves  exist  only  when  the  shear  velocity  in  the 
layer  is  less  than  that  of  the  half-space.  Love  waves  can  exist  in  a  multi-layered  medium  in  which 
velocity  generally  increases  with  depth,  i.e.  even  if  local  low  velocity  layers  are  present  (Pilant 
1979,  p.  161 ;  Ewing  et  al.  1957,  p.  224).  Rayleigh  waves  propagate  without  dispersion  at  the  sur¬ 
face  of  a  half-space  and,  like  Love  waves,  with  dispersion  in  a  layered  medium.  However,  in  the 
case  of  a  layer  (shear  velocity  Pi )  of  thickness  H  over  a  half-space,  the  Rayleigh  waves  are  not  dis¬ 
persed  but  rather  propagate  with  a  velocity  of  Cr  =  0.92(3 1  if  X  «H  for  all  wavelengths  (Press 
and  Ewing  1950). 

The  nature  of  Rayleigh  waves  is  evident  in  the  derivation  of  the  displacement  equations  for 
waves  propagating  along  the  surface  of  an  elastic  half-space.  The  displacement  vector  for  a  propa¬ 
gating  wave  is  u.  A  property  of  vectors  is  that  any  vector  can  be  written  as  the  sum  of  a  gradient 
and  a  curl;  therefore,  u  =  v  +  w  where  v  =  grad  <p  and  w  =  curl  i p.  The  potentials  <p  and  \p  can  be 
expressed  in  general  terms  as* 

<p  =  A  exp  iu>  [/  -(x/c)  +az] 

•F.A.  Okal,  Yale  University,  personal  communication,  1979. 


are  considered;  layers  j  >  m  compose  the  reflecting  zone  for  which  the  complete  complex  reflec¬ 
tivity  is  Rpp  (oj,  7);  u,  =  \Zk2a  -  k*a  sin2  7;  and  range  A  is  the  angular  separation  between  source 
and  receiver  positions  measured  from  the  center  of  the  earth  (horizontal  range  at  the  surface  is 
*  =  0.017  Are  where  re  is  the  radius  of  the  earth. 

In  evaluating  displacement  with  this  method,  a  Bessel  function  typically  is  replaced  by  its 
asymptotic  approximation  for  large  arguments;  this  is  valid  for  large  source-receiver  distances  or 
large  horizontal  wave  numbers,  i.e.  high  frequencies.  The  integration  is  usually  carried  out  over 
a  limited  range  (due  to  computational  costs)  of  real  ray  parameters.  The  use  of  real  ray  param¬ 
eters  restricts  the  method  to  body  waves;  the  use  of  a  limited  range  of  ray  parameters  introduces 
a  truncation  phase  in  the  seismogram .  The  range  of  ray  parameters  corresponds  to  an  equally  re¬ 
stricted  range  of  angles  of  incidence;  in  the  case  of  a  deep  reflecting-refracting  zone,  horizontally 
traveling  energy  is  left  out  of  the  computation.  If  the  top  of  the  reflecting-refracting  zone  is  at 
the  surface,  horizontally  traveling  energy  contained  in  reflected  and  refracted  waves  as  well  as 
surface  wave  energy  should  be  included  in  the  computation. 

Burdick  and  Orcutt  (1979),  in  a  comparison  of  the  reflectivity  and  generalized  ray  methods, 
cite  the  approximations  customarily  made  and  reference  papers  in  which  such  approximations 
are  not  made  in  obtaining  a  solution.  Sato  and  Hirata  (1980)  evaluate  results  obtained  with  these 
methods  and  with  a  third  hybrid  method  having  characteristics  of  both. 

An  approach  similar  to  the  reflectivity  method  is  given  by  Kennett  (1979, 1983)  for  plane 
waves  from  a  surface  source  that  are  multiply  reflected  during  propagation  in  a  multi-layered 
medium;  effects  of  excluding  converted  shear  waves  and  free  surface  reflections  and  of  including 
attenuation  are  shown. 

The  final  wave  theory  considered  is  full  wave  theory.  This  is  actually  a  term  that  applies  to 
several  theories.  The  only  distinction  among  the  theories  given  here  is  that  one  approach  to  a  full 
wave  theory  of  body  wave  propagation  results  in  the  calculation  of  the  WKBJ  seismogram. 

WKBJ  theory  refers  to  a  method  of  finding  approximate  solutions  to  a  second-order  differ¬ 
ential  equation.  When  the  method  is  applied  to  the  wave  equation,  it  results  eventually  in  the 
WKBJ  seismogram.  Several  approximations  are  made  in  the  process  of  obtaining  the  WKBJ  solu¬ 
tion,  one  of  which  is  to  retain  only  the  first  two  terms  of  a  binomial  expansion.  As  a  consequence 
of  the  approximations,  the  solution  is  valid  for  frequencies  that  are  large  in  relation  to  velocity 
gradients  or  equivalently,  for  a  nearly  homogeneous  medium.  Grant  and  West  (1965)  show  the 
equivalence  between  this  asymptotic  solution  to  the  wave  equation  and  the  approximate  solution 
obtained  by  means  of  the  eikonal  equation,  both  of  which  involve  the  propagation  of  wave  fronts 
along  paths  that  are  determined  by  Snell’s  law.  The  WKBJ  method  cannot  handle  the  situation 
of  a  turning  point,  i.e.  when  the  ray  has  reached  its  maximum  penetration  depth  and  is  directed 
horizontally  prior  to  turning  toward  the  surface,  due  to  a  singularity  that  arises  from  the  use  of 
the  truncated  binomial  expansion.  A  different  method,  such  as  an  Airy  function,  is  used  in  the 
vicinity  of  the  turning  point.  The  solutions  obtained  are  oscillatory  above  the  turning  point  and 
exponentially  decaying  below  the  turning  point. 

The  full  wave  theory  (Aki  and  Richards  1980,  Choy  et  al.  1980)  is  applicable  to  body  wave 
propagation  in  situations  involving  radial  inhomogeneity,  smoothly  varying  density  and  velocity, 
ray  paths  with  turning  points  (including  those  near  grazing  incidence  at  an  interface),  critically 
incident  rays,  caustics,  wave  interference  due  to  multipathing,  and  non-ray  theoretical  effects 
such  as  diffraction.  As  a  result  of  the  first  two  characteristics,  a  realistic  earth  model  can  be  used, 
rather  than  a  representation  of  planar  layers  of  constant  velocity  and  density.  The  great  versatility 
of  the  full  wave  theory  results  from  the  incorporation  of  the  frequency  dependence  of  the  waves. 
The  equation  for  radial  displacement  is  integrated  in  the  real  frequency  domain  and  in  the  ray 
parameter  domain;  depending  on  the  path  of  integration  in  the  ray  parameter  domain,  the  dis¬ 
placement  equation  reduces  to  the  solution  obtained  from  ray  theory  if  the  reflection  and  refrac¬ 
tion  coefficients  are  assumed  to  be  independent  of  frequency.  Chapman  (1978)  has  investigated 
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with  amplitudes  lower  than  a  preselected  percentage  of  the  maximum  amplitude  of  the  first  ar¬ 
rivals  can  be  excluded  from  the  seismogram.  Krebes  and  Hron  (1981)  show  that  there  is  good 
agreement  between  synthetic  seismograms  obtained  using  asymptotic  ray  theory  and  those  ob¬ 
tained  using  the  Thomson-Haskell  matrix  method  in  the  case  of  anelastic  media.  The  method  can 
be  used  with  laterally  varying  structures  (McMechan  and  Mooney  1980)  provided  that  the  rays  can 
be  traced  through  the  structure. 

Quantized  ray  theory  (Wiggins  and  Madrid  1974)  and  subsequently  disk  ray  theory  (Wiggins 
1976)  were  developed  as  modifications  of  asymptotic  ray  theory  to  obtain  more  accurate  body 
wave  amplitudes.  Examples  of  synthetic  seismograms  obtained  with  the  disk  ray  method  for 
models  of  laterally  heterogeneous  elastic  and  anelastic  media  are  given  by  Braile  et  al.  (1983). 

The  third  wave  theory  considered  is  the  reflectivity  method  (Fuchs  1968;  Fuchs  and  Muller 
1971;  Aki  and  Richards  1980,  p.  393  404).  It  is  used  to  study  the  reflection  and  refraction  (in¬ 
cluding  critical  refraction)  of  spherical  waves  incident  on  and  propagating  within  a  zone  of  layers 
with  planar  interfaces.  The  method  is  customarily  applied  to  body  wave  problems  but  can  in¬ 
corporate  surface  waves.  It  involves  the  determination  of  the  progressive  partitioning  of  wave 
energy  due  to  one  or  multiple  reflections  of  each  wave  considered  together  with  phase  changes 
(delays)  due  to  the  propagation  of  the  wave  through  each  layer.  The  layered  zone  in  which  re¬ 
flections  and  refractions  of  interest  occur  can  be  located  at  the  surface  (Fertig  and  Muller  1978) 
or  buried;  in  the  latter  case,  waves  from  a  source  position  may  be  propagated  to  the  reflecting- 
refracting  zone  and  returned  from  the  zone  to  the  receiver  position  by  considering  only  the  oc¬ 
currence  of  refraction  (Fuchs  and  Muller  1971,  Aki  and  Richards  1980)  or  reflection  and  conver¬ 
sion  of  waves  propagating  in  the  surface  region  (Kohketsu  1981).  An  equation  for  displacement 
at  the  receiver  position  is  obtained  after  two  integrations,  the  first  in  the  domain  of  the  wave 
number  (k),  ray  parameter  (p  =  k/oS)  or  angle  of  incidence  and  the  second  in  the  frequency  (go) 
domain.  The  reflectivity  term  in  the  displacement  equation  is  the  exact  sum  of  the  infinite  series 
of  reflections  for  all  possible  multiply  reflected  waves  (including  converted  waves)  within  the  re- 
flecting-refracting  zone  (Aki  and  Richards  1980,  p.  402).  Accordingly,  the  solution  obtained 
with  the  reflectivity  method  is  equivalent  to  the  net  displacement  due  to  the  arrival  of  an  infinite 
number  of  generalized  rays. 

The  compressional  potential  obtained  from  the  reflectivity  theory  for  the  case  of  a  point  source 
at  the  surface  of  a  plane  layered  earth  is 


A 

<P  (w,  A)  =  F(oj)  J  kat 


sm7  cos7 


/yi  \/2nMa  sin7 


exp(4(MQmsin7  -  4))] 


•  T( go,  7)  •  Kpp  (w-  7)  exp(-i'r)  dy 

(Burdick  and  Orcutt  1979)  where 

/•'(go)  =  the  Fourier-transformed  source  time  function 
7  =  angle  of  incidence  in  the  mth  layer 

ka,  =  w/tXj  . 

The  index  /  refers  to  the /th  layer;  layers /  =  1,2  ...,  m  compose  the  refracting  zone  in  which  only 
the  transmission  losses  T( go,  7)  and  the  phase  shifts  r(eo,  7)  of  a  plane  wave  at  constant  7  and  w 
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in  order  to  reduce  computational  time;  here  r  is  the  travel  time  of  the  ray  with  ray  parameter  p 
and  the  potential  equation  is 


where  R  =  \Jr2  +  z2  .  Similar  equations  for  transmitted  waves  and  reflected  waves  in  a  layered 
medium  involve  the  potential  transmission  coefficient  and  the  potential  reflection  coefficient, 
respectively,  and  an  R  term  that  is  dependent  on  the  location  of  the  source  relative  to  the  layer 
boundaries. 

The  displacement  field  at  the  receiver  position  is  determined  by  summing  in  the  time  domain 
the  contributions  of  rays  that  travel  from  the  source  to  that  point.  Each  reflection  or  transmission 
of  a  ray  at  a  layer  interface  is  taken  into  account  as  it  proceeds  from  source  to  receiver.  Computa¬ 
tional  cost  may  impose  a  limit  on  the  number  of  rays  included  in  the  construction  of  the  seismo¬ 
gram.  A  first  attempt  is  to  use  only  the  rays  that  are  reflected  once  between  source  and  receiver; 
the  number  of  such  rays  depends  on  the  number  of  layers  used  to  model  the  velocity  profile  of  the 
medium.  As  multiply  reflected  rays  of  increasing  order  are  included  in  the  computations,  the  syn¬ 
thetic  seismogram  changes  until  eventually  inclusion  of  additional  rays  does  not  further  change  the 
character  of  the  seismogram.  For  a  medium  composed  of  thin  layers,  the  number  of  rays  that  must 
be  included  for  accuracy  is  high;  in  this  case,  truncation  of  the  infinite  series  of  multiple  reflections 
can  introduce  serious  errors.  The  multiples  that  must  be  included  are  specific  to  the  problem  and 
depend  on  the  depth  and  thickness  of  the  reflected  layers,  the  velocity  and  density  contrasts  at 
interfaces,  and  the  source-receiver  separation.  The  validity  of  the  synthetic  seismogram  genera¬ 
ted  depends  upon  the  incorporation  of  all  necessary  rays.  The  solution  is  more  difficult  to  obtain 
if  attenuation  is  incorporated  in  the  problem. 

The  second  wave  theory  considered  is  asymptotic  ray  theory  (Hron  and  Kanasewich  1971 , 

Pilant  1979).  This  approach  represents  the  displacement  due  to  propagation  of  a  plane  wave  in  a 
layered  medium  as  the  sum  of  an  infinite  series  of  products  involving  a  spatially  dependent  ampli¬ 
tude  and  two  frequency  terms,  one  an  exponential  and  one  a  power  series  in  inverse  frequency, 


oo  eiu(t-T) 

u  (x,  y,  z,  0  =  E  An(x,  y,  z) -  ; 

"=0  (ioj)n 

here  T  is  a  phase  function  that  depends  upon  the  ray  path.  The  result  is  an  asymptotic  solution 
to  the  wave  equation  that  has  characteristics  of  geometrical  ray  theory.  In  the  limit  of  high  fre¬ 
quencies,  a  useful  asymptotic  solution  is  obtained  with  only  one  or  two  terms  retained.  If  only 
the  first  term  («=0)  of  the  series  is  retained,  the  displacement  equation  is  the  same  as  the  solution 
to  the  eikonal  equation  (eq  4)  and  the  results  obtained  are  necessarily  identical  with  those  due  to 
geometrical  ray  theory  for  wave  propagation  in  a  plane-layered  medium.  Higher-order,  frequency- 
dependent  terms  are  used  in  problems  involving  diffraction  phenomena,  curved  interfaces,  and 
inhomogeneous  media;  the  n=  1  term  is  necessary  in  the  case  of  head  waves. 

Hron  et  al.  (1974)  give  examples  of  synthetic  seismograms  generated  by  the  asymptotic  ray 
method  for  the  case  of  a  compressional  plane  wave  impinging  on  a  plane-layered  medium  and  a 
surface  source.  Each  ray  is  specified  by  the  number  of  linear  segments  composing  it  and  the  num¬ 
ber  of  conversions  that  occur.  An  advantage  of  this  method  is  that  each  constitutent  ray  (or  group 
of  kinematic  rays)  of  the  synthetic  seismogram  can  be  identified.  The  amplitude  of  each  ray  is 
determined  using  plane  wave  coefficients  of  reflection  and  transmission;  the  use  of  frequency- 
independent  coefficients  is  justified  by  the  use  of  zero-order  terms  (high  frequency  approxima¬ 
tion).  Geometrical  spreading,  neglected  here  because  the  source  is  assumed  to  be  at  infinity,  is 
normally  incorporated  in  the  ray  amplitude  formulation  (McMechan  and  Mooney  1980).  Rays 


formed  by  the  intersection  of  rays)  with  finite  amplitudes  rather  than  the  infinite  amplitudes  that 
ray  theory  predicts,  and  the  interference  of  waves  that  travel  different  paths  but  have  similar  ar¬ 
rival  times  (multipathing).  The  requirement  concerning  velocity  changes  (X0  ( 5C'/C)  «  1)  is 
satisfied  it  the  wave  is  of  sufficiently  short  wavelength  (high  frequency  approximation)  or  if  abrupt 
velocity  changes  in  the  actual  velocity  profile  are  replaced  by  layered  transition  zones  of  constant 
velocity  in  which  the  thickness  of  the  layers  is  such  that  the  requirement  is  not  violated,  i.e.  thick¬ 
ness/wavelength  >  10  (or  even  100).  The  requirement  concerning  changes  in  curvature  of  the  wave 
front  (X06  b' "  «  1 ),  which  again  is  satisfied  in  the  limit  of  high  frequencies,  limits  the  validity  of 
the  ray  theory  approach  to  plane  waves;  this  requirement  is  a  restriction  to  the  far-field  region  such 
that  spherical  waves  have  traveled  far  enough  from  the  source  so  as  to  appear  as  plane  waves.  The 
far-field  condition  requires  that  the  horizontal  source-receiver  separation  be  equal  to  several  wave¬ 
lengths. 

Wave  theory 

There  are  several  methods  of  generating  synthetic  seismograms  that  are  based  upon  wave  theory. 
In  each  case,  the  wave  theory  can  he  shown  to  encompass  geometrical  ray  theory.  A  significant 
difference  between  ra>  theory  and  wave  theory,  which  often  also  formulates  its  development  in 
terms  ot  rays,  is  that  wave  theory  incorporates  a  dependence  on  the  frequency  of  the  seismic  wave 
in  its  solution. 

The  first  wave  theoiy  considered  is  generalized  ray  theory,  also  known  as  exact  ray  theory,  the 
method  ot  generalized  reflection  and  transmission  coefficients  (Muller  1970)  or  the  Cagniard-de 
Hoop  method  (Aki  and  Richards,  1 980,  224-253,  386-393).  It  is  applicable  to  problems  involving 
spherical  waves  incident  on  planar  interfaces  and  can  incorporate  direct  waves,  reflected  waves, 
critically  refracted  waves  and  Rayleigh  waves  in  the  solution.  It  requires  an  integration  in  the  com¬ 
plex  ray  parameter  plane  and  a  Laplace  transform  to  obtain  solutions  in  the  time  domain;  the  path 
of  integration  is  different  for  each  generalized  ray  and  must  be  found  exactly.  In  Laplace  space, 
the  cornpressional  potential  due  to  a  point  source  in  a  homogeneous  medium  is+ 
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where  p  =  the  complex  wave  parameter 

s  =  an  imaginary  parameter  analogous  to  the  frequency  of  a  wave  ( p  =  y/s,  where  y  is 
analogous  to  a  radial  wave  number) 
a  =  X  +  2p 
r?a  =  0/«2  -P2)Vl 
G  =  Green’s  function 
A'0  =  the  modified  Bessel  function. 

No  approximations  need  be  made  in  obtaining  the  solution.  If  Bessel’s  function  is  replaced  by 
the  first  term  of  its  asymptotic  expansion,  the  expression  for  0  in  the  time  domain  that  results  is 
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where  H  is  the  Heaviside  function  and  *  denotes  convolution^ .  The  first-motion  approximation 
(t-T  +  Ipr)^1 (2 pr)x/l,  which  is  valid  for  large  ranges  (r)  but  fails  at  long  times  (f),  is  often  made 
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i.  e.  when  sin0  =  1 ;  this  is  the  penetration  depth  of  the  ray  which  then  turns  upward.  The  time  re¬ 
quired  for  the  ray  to  travel  from  depth  Z0  to  depth  Zf  is  given  by 


zf 

,  _ ,  _  f  dZ 
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(6) 


The  horizontal  (surface)  component  of  the  ray  path  at  a  depth  Zf  is 


'  -JLCl&tZ  . 

{  [1  -P2C\Z)]* 

^  A 


(7) 


If  the  ray  reaches  Zpthen  the  integration  of  eq  6  and  7  must  be  performed  twice,  once  for  the 
downgoing  ray  path  and  once  for  the  upgoing  ray  path.  The  accuracy  of  all  the  calculations  de¬ 
pends  both  on  the  approximations  of  ray  theory  stated  above  and  on  the  velocity  function  C(Z) 
used  to  approximate  the  actual  variation  of  velocity  with  depth.  Greenhalgh  and  King  (1981 )  give 
the  ray  equations  that  are  valid  for  various  velocity -depth  functions.  Also,  the  equations  were  de¬ 
rived  for  the  case  of  a  flat  earth  that  is  homogeneous  or  is  inhomogeneous  with  depth  only.  The 
approximation  of  a  flat  earth  is  justified  in  wave  propagation  problems  involving  source-receiver 
separations  small  enough  that  the  curvature  of  the  earth  may  be  neglected.  Officer  (1974)  derives 
equivalent  equations  for  the  case  of  a  spherically  stratified  earth;  for  this  case,  the  quantity  that  is 
constant  along  the  ray  is  r  sin  0/C,  which  is  also  called  the  ray  parameter,  where  r  is  the  radial  dis¬ 
tance  from  the  center  of  the  earth. 

Shadow  zones  exist  when  the  velocity  profile  is  such  that  rays  are  directed  away  from  a  parti¬ 
cular  region;  this  occurs  when  a  trend  of  gradually  increasing  velocity  with  depth  is  interrupted  by 
a  negative  velocity  gradient  (velocity  decreases  with  depth)  and  then  resumed.  Similarly,  conver¬ 
gence  of  rays  occurs  when  the  velocity  profile  is  such  that  rays  cross  and  become  focused  in  a  re¬ 
gion;  in  this  case,  a  trend  of  gradually  increasing  velocity  with  depth  is  interrupted  by  a  greater 
positive  velocity  gradient  and  then  resumed.  A  complete  discussion  of  variation  in  the  location 
and  energy  content  of  rays  arriving  at  the  surface  due  to  different  velocity  profiles  is  given  in 
Officer  (1974). 

The  path  that  a  given  ray  travels  through  a  medium  can  be  determined  by  ray  tracing;  this  re¬ 
quires  specifying  the  takeoff  angle  (the  orientation  at  which  the  ray  leaves  the  source)  and  the 
velocity  profile  through  the  medium  and  then  applying  Snell’s  law  sequentially.  Repeating  this 
process  for  a  range  of  takeoff  angles  leads  to  identification  of  locations  at  which  wave  energy  ar¬ 
rives  at  the  surface,  as  well  as  focusing  and  defocusing  of  the  energy.  The  converse  problem,  de¬ 
termining  what  rays  from  a  known  source  arrive  at  a  given  surface  location,  is  more  difficult  and 
involves  an  iterative  process  beginning  with  a  trial  solution.  In  the  shooting  method,  the  initial 
orientation  of  the  ray  is  specified  and  subsequently  adjusted  until  the  ray  arrives  acceptably  close 
to  the  receiver  location;  in  the  bending  method,  the  source  and  receiver  positions  are  specified  and 
the  ray  path  is  adjusted  (Aki  and  Richards  1980,  p.  727).  Procedures  for  ray  tracing  when  source 
and  receiver  positions  are  known  are  given  by  Chander  (1977);  they  are  applicable  to  the  situa¬ 
tion  of  dipping  planar  layers.  Because  these  methods  yield  frequency-independent  results,  the  in¬ 
terference  of  waves  arriving  at  the  same  time  is  found  by  determining  the  phase  change  along  each 
ray  path  and  summing  the  results. 

Limitations  on  the  usefulness  of  geometrical  ray  theory  result  from  approximations  that  re¬ 
strict  its  validity  to  situations  in  which  changes  in  velocity  or  curvature  are  small  over  a  wave¬ 
length,  and  from  its  inability  to  predict  the  following  frequency-dependent  effects:  the  arrival 
of  body  waves  in  shadow  zones  due  to  diffraction,  the  arrival  of  waves  near  a  caustic  (the  envelope 


9 


where  E  =  the  amount  of  energy  per  unit  time  per  unit  solid  angle  emitted 
by  a  point  source  at  the  origin  of  the  coordinate  system 
60  =  the  angle  with  respect  to  the  z-axis  at  which  the  ray  leaves  the  source 
.v  =  the  .v -component  of  the  distance  along  the  ray  path. 

The  decrease  in  intensity  due  to  geometrical  (cylindrical)  spreading  is  1  fx\  the  change  in  intensity 
due  to  changes  in  the  cross-sectional  area  of  the  ray  bundle  caused  by  focusing  is  tan  80l(dx/d60) 
(Officer  1974,  p.  214).  The  assumption  that  energy  does  not  cross  rays  implies  that  energy  is 
neither  diffracted  according  to  wave  theory  nor  scattered  by  inhomogeneities  along  the  ray.  The 
presence  of  wave  energy  in  shadow  zones,  locations  for  which  ray  theory  predicts  no  surface  ar¬ 
rivals,  is  attributed  to  scatterers  and  the  occurrence  of  diffraction  that  direct  energy  into  the 
shadow  zone  (Urick  1983). 

A  ray  in  a  homogeneous  medium  is  traced  using  Snell’s  law, 
sin0o/C  =  p 

where  p  =  constant  for  a  given  ray  and  is  called  the  ray  parameter 

=  the  angle  the  ray  makes  with  the  vertical  to  the  surface  (takeoff  angle) 

C  =  the  wave  velocity  in  the  medium. 

If  wave  velocity  changes  along  the  ray,  the  angle  9  at  which  the  ray  is  directed  must  change  in 
order  for  p  to  remain  constant.  This  relationship  is  given  as 


for  an  inhomogeneous  medium  modeled  as  n  layers  of  constant  velocity  Cn;  0n  is  the  angle  indi¬ 
cating  the  orientation  of  the  ray  in  a  layer  of  velocity  Cn.  As  the  thickness  of  the  layers  goes  to 
zero,  C  is  a  continuous  function  of  depth,  C(Z),  and  eq  5  can  be  written  as 

si11  _  sin 

c(4)  =  c(Z)  =  p 

(Clay  and  Medwin  1977,  p.  83;  Telford  et  al.  1976). 

There  is  additional  information  on  the  propagation  of  energy  that  can  be  obtained  from  ray 
theory  (Officer  1974,  Clay  and  Medwin  1977).  The  relationship  between  curvature  of  the  ray, 
dO/dS,  and  velocity  gradient,  dQdZ ,  is 

dd_=  dC 
dS  PdZ 

where  5  is  the  coordinate  distance  measured  along  the  ray.  The  ray  is  a  straight  line  in  a  medium 
of  constant  velocity;  it  is  a  portion  of  a  circle  in  a  medium  with  a  linear  velocity  profile,  in  which 
case  the  ray  curves  toward  the  region  of  minimum  velocity.  The  ray  becomes  horizontal  at  a  depth 


Solutions  ot  the  eikonal  equation  have  the  form  W(x,  y,  z)  is  constant  and  represent  wave  fronts. 
At  a  given  time  t,  all  points  along  a  particular  wavefront  W  will  be  in  phase.  A  normal  to  the  wave 
front  along  a  particular  wave  front  defines  the  direction  of  propagation  and  is  a  ray.  Officer 
(1974)  shows  that 

U  -  A(x,  y,  z)  e,tJ  I  w(x,y.z)/c0  - 1 1 


then  W  is  a  solution  to  the  eikonal  equation;  here,  A0  =  2itC0/u>  is  the  wavelength  of  the  distur¬ 
bance  at  angular  frequency  to  in  the  reference  medium  of  velocity  C0.  This  condition  is  satisfied 
in  the  limit  of  high  frequencies  (A->-0).  Of  practical  interest,  W  is  a  solution  to  the  eikonal  equa¬ 
tion  and  U  is  a  good  approximation  to  a  solution  of  the  wave  equation  if  A02  (A" /A) « ( W ')2 
(prime  denotes  differentiation  with  respect  to  a  spatial  coordinate).  This  condition  can  be  re¬ 
written  to  show  that  a  solution  to  the  eikonal  equation  will  satisfactorily  approximate  a  solution 
to  the  wave  equation  if  1)  the  fractional  change  in  velocity  gradient,  6 C',  over  a  wavelength  is 
small  compared  with  the  gradient  Cl\,(\0  ( dC'/C )  «  1),  2)  the  change  of  curvature  of  the 
wave  front  is  small  over  a  wavelength,  (A06W"  «  1),  or  3)  the  fractional  change  in  the  space  rate 
of  change  of  amplitude  is  small  over  a  wavelength  (\q  (6A  ’/A)  «  1)  (Officer  1974,  p.  205). 

The  direction  in  which  wave  energy  propagates  is  parallel  to  the  rays.  Because  it  is  assumed 
that  energy  does  not  cross  rays  (i.e.  the  energy  contained  within  a  bundle  of  rays  [raytube]  at  the 
source  remains  within  the  space  defined  by  those  rays)  the  propagation  of  energy  through  a  medi¬ 
um  becomes  known  by  determining  the  paths  of  rays  through  the  medium  (Fig.  6).  The  space  de¬ 
fined  by  a  given  bundle  of  rays  changes  due  to  geometrical  spreading  and  due  to  changes  in  veloc¬ 
ity  with  depth  that  result  in  focusing  or  defocusing  of  the  rays.  Accordingly,  the  intensity  of  the 
disturbance  varies  along  the  ray  paths.  The  intensity  at  a  point  along  a  ray  path  between  source 
and  receiver  at  the  same  depth  is  given  by 

E  tan0Q 
/(jC,6,o)  =  x(dx/dd0) 


Figure  6.  Rays  A  and  B  emanating  from  a 
source  Sat  takeoff  angles  of  6  0  and  0o+d0o, 
respectively.  The  rays  are  linear  in  a  medium 
of  constant  velocity,  V.  In  time  At,  each  ray 
travels  a  distance  D  =  VAt.  The  location  of 
ray  A  and  the  wave  energ}’  associated  with 
it  is  then  (X0  +  D  sinOQ,  Z0  +  D  cos60J. 


When  the  angle  of  incidence  of  a  wave 
at  an  interface  is  real,  the  coefficients  (trans¬ 
mission,  reflection)  for  plane  waves  at  planar 
interfaces  are  real  and  independent  of  fre¬ 
quency.  The  latter  property  means  that 
there  is  no  change  in  wave  shape  at  the  in¬ 
terface;  i.e.  the  departing  transmitted  and 
reflected  waves  have  the  same  shape  as  the 
incident  wave  but  probably  different  ampli¬ 
tudes.  If  either  the  incident  wave  or  the  in¬ 
terface  is  nonplanar,  there  may  be  a  change 
in  pulse  shape  (a  phase  change)  at  the  inter¬ 
face  (Aki  and  Richards  1980,  p.  1 55). 

The  angle  at  which  plane  waves  are  inci¬ 
dent  at  a  boundary  is  crucial  to  the  long- 
range  propagation  of  wave  energy  within  a 
layer  or  a  waveguide  (Fig.  5).  When  the 
angle  of  incidence  is  less  than  the  relevant 
critical  angle,  wave  energy  is  lost  across 
each  boundary  in  the  form  of  transmitted 


Figure  5.  Reflected  waves  and  transmitted 
waves  generated  by  plane  waves  incident  at 
the  waveguide  boundaries.  At  angles  less 
than  critical  these  are  represented  by  solid 
rays.  Reflected  waves  and  head  waves  gener¬ 
ated  by  plane  waves  that  are  critically  inci¬ 
dent  at  the  waveguide  boundaries  are  repre¬ 
sented  by  dashed  rays.  Reflected  waves  gen¬ 
erated  by  plane  waves  incident  at  angles 
greater  than  critical  are  represented  by  dot¬ 
ted  rays;  the  associated  inhomogeneous 
waves  are  represented  by  wavy  lines. 


( P ,  5)  waves.  When  the  waves  are  incident 

at  the  critical  angle,  reflected  waves  and  a  head  wave  exist  and  the  wave  energy  is  effectively  trapped 
by  the  boundaries  of  the  waveguide.  When  the  angle  of  incidence  is  greater  than  the  critical  angle, 
the  criterion  for  total  reflection  is  satisfied  but  some  energy  is  lost  across  the  boundary  from  the 
inhomogeneous  waves.  The  loss  of  energy  from  the  waveguide  is  significantly  lower  when  waves 
are  incident  at  the  waveguide  boundaries  at  angles  greater  than  critical. 

The  wave  energy  in  a  layer  or  waveguide  may  be  described  by  a  sum  of  normal  modes.  Each  nor¬ 
mal  mode  represents  horizontally  propagating  waves  that  have  the  profile  of  standing  waves  in  the 
z-direction;  it  is  a  solution  to  the  wave  equation  and  satisfies  the  boundary  conditions  of  the  prob¬ 
lem.  The  stress-depth  profile  of  a  normal  mode  may  be  determined  from  an  analysis  of  the  con¬ 
structive  interference  of  plane  waves  that  are  totally  reflected  within  the  waveguide. 


SECTION  3.  BODY  WAVES:  RAY  THEORY  AND  WAVE  THEORY 

Ray  theory  and  wave  theory  are  each  the  basis  of  analytical  approaches  to  the  study  of  the 
propagation  of  body  waves  in  the  earth.  Their  similarity  is  that  each  yields  asymptotic  solutions 
to  the  wave  equation;  their  differences  lie  in  the  approximations  used  in  their  development  and, 
consequently,  in  the  restrictions  on  their  validity. 

Geometric  ray  theory 

The  concept  of  rays  arises  from  the  eikonal  equation, 


(irS  +  dr)2  +  (I? )'  =%t="2’ 


dx'  vdy 


where  C  =  velocity 

C0  =  a  constant,  reference  velocity 
n  =  the  index  of  refraction  (C0/C  =  n). 
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where  A0  is  an  arbitrary  constant  used  to  define  the  wave  number  A 2  =  n2kQ  ,  n  being  the  index  of 
refraction.  If/ is  assumed  to  vary  slowly  with  respect  to  range,  then  \d2f/dr 2 1 «  |2 k0  (9//3r)|  and 
the  Helmoltz  equation  is  simplified  to  the  parabolic  wave  equation  (d2f/dz2 )  +  2ik0(df/dr)  +  (A2  - 
A'02)/=  0.  Computer  codes  that  solve  the  parabolic  equation  also  are  used  in  studying  problems  in 
underwater  acoustics.  The  solution  to  the  Helmholtz  or  the  parabolic  equation  for  particular  bound¬ 
ary  conditions  (sound  speed  profile,  source  depth,  receiver  depth)  is  presented  as  a  plot  of  propaga¬ 
tion  loss  as  a  function  of  range  for  a  specified  frequency.  Coppens  ( 1 982)  discusses  the  use  of  the 
parabolic  equation  and  cites  errors  in  the  calculated  acoustic  field  that  result  from  the  assumption 
invoked  in  deriving  it.  One  error  is  that  the  phase  velocity  of  a  particular  normal  mode  obtained 
by  solving  the  Helmholtz  equation  is  different  from  the  equivalent  term  obtained  by  solving  the 
parabolic  equation,  and  so  the  spatial  pattern  of  the  phase-coherent  combination  of  the  pressure 
terms  will  be  different  for  the  two  solutions.  The  parabolic  equation  may  be  used  under  the  con¬ 
ditions  that  the  fractional  change  in  velocity  over  one  wavelength  is  small,  and  the  angle  of  propaga¬ 
tion  with  respect  to  the  horizontal  is  small  (Urick  1979). 

When  an  integration  over  wave  number  is  done  to  obtain  an  integral  solution  to  the  wave  equa¬ 
tion,  a  contour  of  integration  along  the  real  A-axis  encounters  no  poles  due  to  normal  modes  if  the 
medium  is  assumed  to  be  lossy.  Tire  Fast  Field  Program  (FFP)  uses  the  Fast  Fourier  Transform  to 
numerically  integrate  the  solution  to  the  wave  equation  for  the  situation  of  lossy  layers  over  a  half¬ 
space  and  then  removes  the  effect  of  attenuation  due  to  anelasticity.  The  result  is  values  of  pres¬ 
sure,  particle  displacement,  or  velocity  at  range  points  rn  =  rQ  +  nAr,  n  =  0,  1 ,  ...,  N-l  .  Which 
modes  are  included  in  the  solution  (normal,  leaky)  is  determined  by  the  range  of  wave  numbers 
Am  =  A0  +  m AA,  m  =  0,  l , ...,  yV-1  selected.  Kutschale  (1 973)  notes  that  use  of  the  FFP  is  quicker 
than  use  of  normal  mode  theory  in  obtaining  a  solution  to  the  wave  equation  because  FFP  solves 
the  integral  solution  to  the  wave  equation  directly  while  normal  mode  theory  requires  that  the 
roots  of  the  dispersion  equation  be  determined  and  then  the  corresponding  normal  modes  be 
summed.  Kutschale  uses  the  FFP  method  to  investigate  how  several  factors  (presence  of  a  surface 
ice  layer,  rigidity  of  bottom  sediments,  source  depth,  small  changes  in  sound  speed  profile  in  upper 
400  m  of  water,  variations  in  ice  roughness)  affect  compressional  wave  propagation  in  the  ocean. 


SECTION  6.  FINITE-DIFFERENCE  METHOD 

The  basis  of  the  finite-difference  method  is  that  it  provides  an  exact  solution  to  an  approxima¬ 
tion  of  a  differential  equation.  In  seismology  it  is  the  wave  equation  that  is  approximated.  Where¬ 
as  an  analytic  method  provides  a  continuous  solution  to  the  wave  equation,  u  =  u{x,  z,  t),  a  numer¬ 
ical  method  provides  a  discrete  solution,  u  =  u(mAx,  nAz,pAt),  to  an  approximation  of  the  wave 
equation.  The  domain  of  the  finite -difference  solution  is  shown  schematically  in  Figure  12  as  a 
distribution  of  points  in  which  thp  spacing  between  adjacent  points  is  Ax  parallel  to  the  x-axis, 

Az  parallel  to  the  z-axis,  and  At  parallel  to  the  r-axis.  An  introduction  to  the  finite-difference 
method  and  its  applications  is  given  in  many  texts  (e.g.  Richtmeyer  and  Morton  1967,  Lapidus 
and  Pinder  1982). 

Boore  (1972)  has  reviewed  the  use  of  the  finite-difference  method  in  the  study  of  seismic  wave 
propagation.  He  cites  the  advantages  of  using  this  method  as  being  the  nature  of  the  computations 
(a  program  to  use  the  finite-difference  method  is  easily  written,  is  adapted  readily  to  a  variety  of 
problems,  and  uses  input  data  that  are  easily  prepared)  and  the  use  of  transient  signals  (because  of 
which  one  computer  run  yields  information  on  many  frequencies).  He  also  cites  advantages  in  the 
presentation  of  results:  displacement  at  a  given  position  as  a  function  of  time  (i.e  the  time  series 
recorded  by  a  seismometer),  or  displacement  at  a  given  time  as  a  function  of  position  (i.e.  the  total 
wave  field).  The  finite-difference  method  is  useful  in  solving  problems  for  which  there  is  no  a;,My- 
tical  solution  or  for  which  the  analytical  solution  is  more  costly  or  less  adaptable  due  to  the  com- 


20 


a.  Finite-difference  scheme  for  problems  of  one-dimensional 
wave  propagation. 


b.  Finite-difference  scheme  for  problems  of  two-dimensional 
wave  propagation. 

Figure  12.  Grids  and  templates  for  finite-difference  method. 

(adapted  from  Lapidus  and  Finder  1982). 

plexity  of  the  physical  model.  Use  of  the  finite-difference  method  is  particularly  advantageous 
with  problems  that  concentrate  on  the  near-field  region  of  sources  or  that  involve  wave  propaga¬ 
tion  through  layered  media  in  the  case  where  the  thickness  of  the  layers  is  on  the  order  of  the 
seismic  wavelengths.  Another  discussion  of  the  finite-difference  method  as  applied  to  seismology 
that  includes  more  recent  references  is  given  by  Aki  and  Richards  (1980).  Examples  of  the  ways 
in  which  the  numerical  solution  to  the  wave  equation  may  be  displayed  are  given  in  Figure  13. 

There  are  two  types  of  finite-difference  schemes,  both  of  which  are  recursive.  To  solve  for  the 
future  displacement  at  a  specified  position,  the  explicit  approach  uses  supplied  or  previously  cal¬ 
culated  values  of  displacement  at  that  position  and  neighboring  positions  for  past  and  present 
times.  The  implicit  approach  again  uses  known  displacements  at  the  specified  position  and 
neighboring  positions,  but  for  past,  present  and  future  times.  In  the  notation  used  here,  the  ex¬ 
plicit  approach  uses  displacements  at  time  levels  (p-1)  and  (p)  to  calculate  displacements  at  time 
level  (p+ 1 ),  while  the  implicit  approach  uses  displacements  at  time  levels  (p-1 ),  (p ),  and  (p+1 )  to 
calculate  displacements  at  time  level  (p+1)  (See  Fig.  12).  Emerman  et  al.  (1982)  investigated  the 
use  of  an  implicit  finite -difference  approximation  of  the  wave  equation  and  found  it  to  be  less 
accurate  than  the  explicit  formulation  for  the  same  size  increments  of  space  and  time.  Braile  et 
al.  (1 983),  however,  compared  four  implicit  formulations  (including  that  of  Emerman  et  al.)  of 
the  scalar  wave  equation  with  the  explicit  formulation  of  Alford  et  al.  (1974);  their  results  in¬ 
dicate  that  one  of  the  implicit  formulations  provides  reasonable  accuracy  with  a  grid  spacing  that 
is  approximately  twice  as  coarse  (4  to  5  grid  points  per  wavelength)  as  that  required  by  the  ex¬ 
plicit  and  other  implicit  formulations,  and  thus  results  in  a  factor  of  four  savings  in  computer 
storage  and  a  factor  of  eight  savings  in  computer  time. 

An  understanding  of  the  rather  complex  finite-difference  equations  given  below  can  be  gained 
by  following  the  derivation  of  finite  difference  equations  in  the  much  simpler  one-dimensional 
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a.  Geometry  and  media  labels  for  four 
materials  joining  along  a  vertical  interface. 
If  the  two  conditions  at  the  bottom  are 
satisfied,  an  analytic  solution  to  the  prob¬ 
lem  of  an  incident  Love  wave  can  be  found. 
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b.  Surface  displacements  at  different  instants 
of  time  for  the  problem  indicated  in  a.  The 
location  of  the  vertical  interface  is  given  by 
the  dashed  line.  The  various  elastic  and  geo¬ 
metric  parameters,  Pi  =3.85,  Pi  =3.0,  Pi  = 
4.75,  p2  =  5.65,  P3  =  4.23,  p3  =  7.44,  p4  = 
5.53,  p4=  8.07,  H  =  35.0,  were  chosen  such 
that  a  large  reflected  wave  would  be  present. 
The  units  of  p  and  p  are  km/s  and  g/cm3. 
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c.  Computer-generated  contour  plots  of  the  spatial  distribution  of  Love-wave  displacements  at 
different  instants  of  time  for  the  problem  illustrated  in  a.  The  results  here  are  derived  from  the 
same  experiment  as  those  in  b.  The  top  edge  of  each  plot  corresponds  to  the  free  surface.  The 
vertical  line  denotes  the  vertical  interface.  The  horizontal  interface  is  shown  in  plot  for  T  =  112. 
Vertical  variation  is  distorted  due  to  variable  grid  spacing  with  depth.  Each  plot  can  be  viewed  as 
a  contour  map  of  the  displacements  at  each  grid  point  for  a  certain  instant  of  time,  in  which  the 
area  between  every  other  contour  is  indicated  by  a  given  symbol. 

Figure  13.  Application  of  the  finite-difference  method  (from  Boore  1970). 


case.  The  Jt-axis  is  divided  into  segments  of  length  AX  that  are  bounded  by  grid  points.  The  grid 
points  are  specified  by  m  AX,  where  m  =  0, 1 , 2, ...,  and  displacement  at  grid  points  is  given  by 
U(mAX)  =  Um. 

Finite-difference  approximations  to  the  partial  derivatives  bU/bX  =  Ux,b2  U/bX 2  =  UXK  are  ob¬ 
tained  from  a  series  expansion  of  U( X).  Two  Taylor  series  for  U(X)  about  the  point  (mAX)  = 
are 
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Equations  14  and  15  are  solved  individually  for  Ux  I  to  obtain 
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Two  approximations  to  Ux  at  Xm  are  obtained  from  eq  16  and  17  by  omitting  terms  in  which  AX 
is  of  the  order  1  or  higher: 
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(18) 
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Equation  18  is  termed  a  forward  difference;  Ux  |  is  determined  by  the  displacements  at  Xm  and 
the  adjacent,  higher  index  grid  point,  Xm+i .  Equation  19  is  termed  a  backward  difference;  Ux  | 
is  determined  by  the  displacements  at  Xm  and  the  adjacent,  lower  index  grid  point,  Xm_j . 

A  central  difference  approximation  to  Ux  I  is  obtained  by  adding  eq  16  and  17  and  solving  for 
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By  subtracting  eq  1 7  from  eq  16  and  solving  for  Uxx  i  ,  the  following  finite-difference  approxima¬ 
tion  to  b2  U/bX 2  is  obtained:  m 
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The  first  term  of  the  series  remainder  is  -((AX)2  / 1 2)  Uxxxx  |  where  Xm  -  AX  <  e  <  Xm  +  AX. 

The  equation  which  Boore  (1972)  used  in  his  study  of  Love  surface  waves  and  SH  body  waves 
is 
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where  is  notation  for  displacement  at  a  point  (m  Ax,  «Az)  at  a  time  pAi  and  (32  =  (pip)2 . 

This  equation  is  a  finite-difference  approximation  to  the  partial  differential  equation  for  displace¬ 
ment  in  the y -direction  in  a  homogeneous  medium,  9 2ujdt2  =  |32(92u/9x2+92u/9z2).  It  is  correct 
to  second  order  in  the  increments.  In  the  case  of  an  inhomogeneous  medium, 
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(Aki  and  Richards  1980,  p.  781). 

The  wave  equation  for  the  situation  of  in-plane  waves  (P-SV,  Rayleigh)  propagating  in  a  homo¬ 
geneous  isotropic  body  is 
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In  vector  form  the  wave  equation  is 
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Aki  and  Richards  (1980)  (citing  Alterman  and  Rotenburg  1969)  give  the  following  finite-difference 
approximation  to  the  wave  equation: 


u(x,  z,  t  +  A/)  =  2u(x  ,  z,  r)-  u(x,  z,  t  -  Af)  +  A  (-^)J  [u(x  +  Ax),  z,  f)  -  2u(x,  z,  f) 
+  u(x  -  Ax.z.f)]  +  [u(x+Ax,z+Az,  t)  -  u(x+Ax,z-Az,t) 

-u(x-Ax,z+Az,/)  +  u(x-Ax,z-Az,0]  +c(^f  [u(x,z+Az,r) - 
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-  2u(.v,  z,  t)  +  u(x,  z  -  Az,  f )] . 

This  finite-difference  equation  is  correct  to  second  order  in  the  increments.  Aki  and  Richards  also 
give  Aiterrnan  and  Karal’s  (1968)  finite-difference  equation  for  spherical  waves  propagating  from 
an  explosive  source  in  a  layered  medium,  and  a  finite  difference  equation  for  one-dimensional 
wave  propagation  in  the  case  where  wide-angle  scattering,  back-scattering,  and  reflection  are  negli¬ 
gible;  the  latter  equation  is  suited  to  the  study  of  wave  propagation  in  a  laterally  varying  medium 
from  a  teleseismic  source. 

The  wave  equation  (constitutive  equation)  can  be  approximated  by  several  different  finite  dif¬ 
ference  equations  that  may  all  approach  the  differential  equation  in  the  limit  that  Ax  (and  Az) 
and  A /  approach  zero,  but  that  may  have  different  solutions  when  the  increments  are  finite.  In 
general,  the  choice  of  a  finite  difference  equation  to  approximate  a  differential  equation  affects 
the  accuracy  and  stability  of  the  numerical  solution.  Equations  18,19  and  20  are  all  finite-differ¬ 
ence  approximations  to  Ux  I  .A  measure  of  the  error  that  results  from  the  finite-difference 

1  m 

approximation  is  given  by  the  first  term  of  the  truncated  series  expansion  of  U(X): 

+  AT  «*m  +  AX(eq  16) 

2  He  '  Xm  -AX<e<Xm  (eq  17) 

-  (^p-2  Uxxx  |  ,  Xm  -  AX  <  e  Xm  +  AX  (eq  20) 

Lapidus  and  Pinder  ( 1982)  list  finite-difference  approximations  to  derivitives  for  the  cases  of  one 
and  two  dependent  variables,  and  list  the  order  of  error  resulting  from  truncation  of  the  series. 

Aki  and  Richards  (1980)  evaluate  the  stability  (i.e.  does  the  error  due  to  approximating  derivitives 
grow  as  the  numerical  solution  proceeds,  or  does  it  remain  bounded?)  of  various  finite  difference 
schemes  for  wave  propagation  in  one  direction,  as  given  by 

9 u  _  J_  9T 
9/  p(x)  9x 

where  u  =  du/dr 

T  =  i.(x)  du/dx 
H(x)  =  X(x)  +  2p(x), /’-waves 
l'(x)  =  p(x),  5- waves. 

They  find  that  a  scheme  using  a  central  difference  for  the  x-derivitive  and  a  forward  difference 
for  the  r-derivitive  is  unstable,  while  a  scheme  using  central  differences  for  both  derivitives  is 
stable  if  restrictions  on  At  and  Ax  are  met.  The  optimum  scheme  that  Aki  and  Richards  consid¬ 
er  involves  staggered  grids,  with  the  grid  for  u  being  displaced  from  the  grid  for  T  by  half  a  grid 
spacing  in  both  x  and  t\  the  use  of  central -difference  approximations  leads  to  an  error  equal  to 
'/■>  of  the  error  of  the  latter  scheme  above  (because  the  error  is  proportional  to  the  square  of  the 
sampling  interval)  while  having  the  same  stability  criterion. 

A  criterion  for  stability  has  the  form  At  <  A x/c,  where  c  is  the  local  wave  velocity.  For  a 
solution  to  be  stable,  the  time  interval  between  adjacent  grid  points  must  be  less  than  or  equal 
to  the  time  required  for  a  wave  to  propagate  over  the  space  interval  between  adjacent  grid  points. 
When  this  criterion  is  met,  w(w+l,  w)  is  uniquely  determined  by  u(m,  n-1 ),  u(m,n+l  )and  u(m- 1 , 
n).  The  stability  condition  for  SH- wave  or  Love  wave  propagation  in  a  uniform  medium  of  n 
space  dimensions  is0A//Ax  <  \l\fn  (Boore  1972).  The  stability  criterion  for  the  two-dimension 
in-plane  (/’-Sfand  Rayleigh)  wave  problem  described  above  is  A t/Ax  <  1  /(or2  +  p2  )Vl  (Aki  and 
Richards  1980,  p.  782). 
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Since  the  finite -difference  method  results  in  a  solution  to  the  constitutive  (wave)  equation  only 
at  net  points,  it  would  seem  logical  to  use  a  net  of  small  grid  size  so  that  die  problem  is  well  re¬ 
solved  spatially  and  temporally.  The  cost  to  obtain  the  numerical  solution,  however,  is  dependent 
primarily  on  the  number  of  calculations  that  must  be  performed  and  on  tire  storage  space  required, 
and  so  is  ultimately  determined  by  the  grid  spacing.  A  compromise  approach  is  to  use  a  fine  grid 
size  early  in  the  numerical  solution,  such  as  near  a  seismic  source,  and  then  to  rezone  the  grid  to  a 
larger  spacing  as  the  problem  progresses,  i.e.  at  net  points  corresponding  to  later  times  and  distant 
positions  (Cooper  1971 ),  A  difficulty  of  the  rezoning  technique  is  that  intensive  (mass  independ¬ 
ent )  quantities  such  as  pressure  or  velocity  that  are  initially  specified  for  the  case  of  the  individual, 
small  grid  zones  must  now  be  appropriately  determined  for  the  larger,  composite  zone.  In  a  study 
of  Love-wave  propagation  in  a  plane-lavered  medium,  Boore  ( 1 972 )  increased  the  grid  spacing 
with  depth  in  order  to  delay  the  arrival  of  waves  reflected  from  the  boundary  of  the  region  of 
interest;  this  modification,  based  on  the  decay  with  depth  of  Love  wave  arrivals,  has  had  mixed 
success.  The  maximum  storage  space  required  by  an  explicit  finite  difference  scheme  is  2(A'.Y) 
(/VZ),  where  NX  and  A'Z  are  the  number  of  grid  points  in  the  x  and  :  directions  (Boore  1972). 
Although  each  time  level  requires  (A’.V)  (.VZ)  storage  spaces  and  although  displacements  at  two 
time  levels  ( p ,  p- 1 )  are  needed  to  calculate  the  displacements  at  the  next  time  level  (p+1 ),  once 
is  calculated  u^'n  is  no  longer  needed  so  storage  locations  for  the  (p- 1 )  data  may  be  used 
for  the  (p+l )  data. 

Application  of  the  finite  difference  method  to  determine  future  displacement  requires  that 
either  the  displacement  at  each  net  point  be  known  at  times  /  =  0  and  i  =  A/,  or  that  velocity  and 
displacement  at  time  t  =  0  be  known.  In  studying  wave  propagation  in  an  inhomogeneous  medium, 
it  is  customary  to  situate  the  seismic  source  in  an  adjacent  laterally  homogeneous  medium  so  that 
displacement  at  t  =  0,  A t  (or  displacement  and  particle  velocity  at  t  =  0)  can  be  determined  from 
an  analytic  solution.  Displacement  is  zero  throughout  the  inhomogeneous  region  until  sufficient 
time  has  passed  (determined  by  distance  from  source  and  velocity  of  homogeneous  region)  for  the 
disturbance  to  arrive  at  the  boundary  between  the  two  regions;  changes  in  displacement  with  time 
are  then  determined  numerically  as  the  disturbance  propagates  through  the  inhomogeneous  region. 
The  theoretical  displacements  at  times  t  =  0  and  r  =  A/  that  are  consistent  with  a  given  surface 
wave  shape  are  synthesized  using  a  Fast  Fourier  Transform  method  (Boore  1970,  Fuyuki  and 
Matsumoto  1980). 

Wave  propagation  problems  studied  by  seismologists  customarily  involve  physical  boundaries 
such  as  a  free  surface,  interfaces  of  layered  media  or  the  perimeter  an  inhomogeneous  region. 
Alterman  and  Loewenthal  ( 1972)  describe  a  technique  for  using  finite-difference  equations  for 
boundary  conditions  in  order  to  obtain  displacements  at  a  free  surface  or  at  an  interface.  In  both 
situations  the  finite  difference  net  is  extended  past  the  physical  boundary  (free  surface,  interface) 
so  that  the  boundary  no  longer  terminates  the  net  superimposed  on  the  half-space  or  layer.  Dis¬ 
placements  at  net  points  on  the  fictitious  line  (the  extension  of  the  net)  are  used  in  the  finite  dif¬ 
ference  equations  of  motion  to  calculate  displacements  up  to  and  including  the  physical  boundary. 
This  approach  is  used  also  by  Kelly  et  al.  ( 1976)  and  Fuyuki  and  Matsumoto  (1980). 

Additional  boundaries  are  introduced  when  the  finite  difference  method  is  used  in  order  to  re¬ 
strict  the  size  of  the  domain  for  which  the  numerical  solution  is  obtained  in  accordance  with  the 
amount  of  computer  storage  available.  These  artificial  boundaries  are  sources  of  reflected  waves, 
which  will  affect  the  numerical  solution  adversely  only  when  sufficient  time  has  passed  for  the  re¬ 
flected  waves  to  have  propagated  into  the  region  of  interest.  Boore  ( 1972)  recommends  that  num¬ 
erical  experiments  using  different  distances  to  the  imposed  boundaries  be  done  in  order  to  define 
the  space-time  region  that  is  free  of  contamination  by  waves  introduced  by  the  presence  of  arti¬ 
ficial  boundaries.  The  artificial  boundaries  can  be  characterized  so  that  wave  energy  is  absorbed 
(Fuyuki  and  Matsumoto  1980)  or  reflections  are  canceled  due  to  the  superposition  of  solutions 
of  opposite  sign  (Smith  1974). 
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Complications  with  a  numerical  solution  that  arise  due  to  the  use  of  a  discrete  grid  are  aliasing 
(false  representation  of  frequency  content)  of  the  continuous  time  series  and  dispersion  of  waves 
propagating  on  a  lattice.  These  complications  can  be  avoided  by  restricting  the  frequency -depend¬ 
ent  information  that  is  extracted  from  the  numerical  solution  (Boore  1972).  The  problem  of 
aliasing  can  be  alleviated  by  regulation  of  the  source  so  that  little  or  no  power  is  present  at  values 
of  wave  number  above  the  Nyquist  wave  number,  which  is  inversely  proportional  to  grid  spacing. 
The  effects  of  lattice-induced  dispersion  can  be  eliminated  by  ignoring  information  at  wavelengths 
for  which  there  are  10  or  fewer  net  points  per  wavelength.  Alford  et  al.  (1974)  showed  that,  with 
an  explicit  second-order  difference  scheme,  there  is  little  dependence  of  phase  velocity  and  group 
velocity  on  stability  ratio  ( cAt/Ax )  for  Av/X  <  0.1 ,  where  X  is  wavelength,  at  the  upper  half-power 
frequency  of  the  source. 

Table  1  lists  examples  of  problems  for  which  synthetic  seismograms  were  generated  using  the 
finite-difference  method. 


Table  1 .  Applications  of  the  finite-difference  method  to  seismology. 


Problems  studied  with  finite- 
difference  method 

Wave  type  (source 

Reference 

Lateral  heterogeneity: 

uniform  layer  of  nonconstant 
thickness  (w/planar  free  sur¬ 
face)  covering  a  homogeneous 
half-space 

Love  wave 

Boore  (1970,  1972) 

Hair  of  layers  over  quarter-space 
with  vertical  interface 

Love  wave 

Boore  (1970,  1972) 

Sloping  layer 

Love  wave 

Boore  (1970,  1972) 

Sloping  interface: 

1)  crust/mantle  interface 

2)  basin 

SH- waves  at  vertical 
incidence 

Boore  et  al.  (1971), 
Boore  (1972) 

Semi-infinite  medium  with  sur¬ 
face  topography 

SH- waves  at  vertical 
incidence 

Boore  (1972) 

90°  wedge  (quarter-space)  in 
otherwise  infinite  homo¬ 
geneous  medium 

Line  source  tangent  to 
apex  of  wedge.  In¬ 
cludes  case  of  dif¬ 
fracted  waves 

Alford  et  al.  ( 1 974) 

Problems  of  interest  in  exploration 
tion  seismology: 
homogeneous  layers  separated 
by  plane,  horizontal  interface, 
includes  case  of  a  weathered 
layer 

quarter-space  embedded  in  a 
half-space 

fluid  (brine,  gas)  saturated, 
100-ft-thick  sandstone  layer 
embedded  in  shale 

Also  heterogeneous  case:  elastic 
constants  specified  at  each 
net  point 

Compressional  (P-)  wave 
line  source  located 
within  upper  layer 

Waves  evident : 
direct  P-waves 
surface  reflected 

P-  waves 
head  waves 

P-waves  converted 
into  5-waves  by 
reflection 

Rayleigh  waves 
diffracted  waves 

Kelly  et  al.  ( 1 976) 

Scattering  at  a  trench 

Rayleigh  wave 

1  uyuki and  Matsumoto 
(1980) 

Layered  half-space:  cylindrical 
coordinates 

Impulsive  line  source  located  in 
interior  of  quarter-space  or  in 
interior  of  three-quarter-space 

Pressure  pulse 

P-waves 

Alterman  and  Loewen- 
thal  (1972) 
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Table  1  (cont’d).  Applications  of  the  finite-difference  method  to  seismology. 

Problems  studied  with  finite- 

_ d  if  fere  nee  m  e  thod _  l\'a  ve  type  I  source  Reference 

Heterogeneous  medium:  spherical  Impulsive  point  source 

coordinates 

Shearing  stress  acting  on  a  crack 
in  a  plane  (deep  focus  earth¬ 
quake) 

Shearing  stress  acting  on  a  crack 
in  a  half-plane  (dip-slip  on  a 
vertical  fault  which  breaks 
free  surface) 

Impulsive  line  source  located  inside 
one  of  two  welded  quarter-spaces 

Dislocation  source  in  a  homogeneous  S- waves,  P-waves/  Braile  et  al.  (1983) 

half-space  Dislocation  source 

Dislocation  source  in  a  laterally 
heterogeneous  velocity  structure 
model;  dislocation  source  is  lo¬ 
cated  on  a  subducting  slab 


Layer  over  half-space 

Love  waves 

Kelly  (1  983) 

Thickening  channel 

Conversion  to  body  waves 

Note:  data  band-pass 

Thinning  channel 

evident 

filtered  to  isolate 

Irregular  layer/half-space 

frequency  compon¬ 

interface 

ents  so  as  to  display 

Channel  w  ith  low  velocity  fill 

modes 

over  subsided  layer  and  sub¬ 


sided  half-space 

Channel  filled  with  same  material 
as  subsided  layer 
Channel  with  high  velocity  fill 
over  subsided  layer  and  sub¬ 
sided  half-space 
Layer  over  two  quarter-spaces 
Layer  over  angular  unconformity 
Periodic  channels  in  upper  portion 
of  layer 


SECTION  7.  FINITE-ELEMENT  METHOD 

There  are  a  number  of  similarities  between  the  finite-element  and  the  finite-difference  methods. 
Both  involve  the  division  of  a  continuous,  two-dimensional  structure  into  discrete  portions  defined 
by  nodes  (or  net  points),  both  yield  numerical  solutions  whose  accuracy  is  dependent  on  the  spac¬ 
ing  between  nodes  relative  to  wavelength,  and  both  assume  that  the  displacements  of  nodal  points 
define  the  complete  displacement  field  for  the  structure.  The  distinction  between  the  methods  is 
that  the  finite-difference  solution  is  obtained  at  discrete  points  by  approximating  the  differential 
equations  that  describe  a  problem  involving  an  elastic  continuum,  while  the  finite-element  solu¬ 
tion  is  obtained  by  approximating  the  elastic  continuum  as  a  composite  of  structural  elements 
whose  internal  displacements  are  known  by  interpolation.  Information  on  the  finite-element 
method  is  found  in  books  on  solving  partial  differential  equations  (e.g.  Lapidus  and  Pinder  1982) 
and  in  examples  of  the  application  of  the  method,  e.g.  structural  mechanics  (Clough  1965)  and 
propagation  of  seismic  waves  (Lysmer  and  Drake  1972.  Smith  1975). 

The  finite  element-method  requires  that  a  structure  be  divided  into  a  finite  number  of  plane 
elements  that  meet  at  nodes.  The  commonly  used  elements  are  quadrilaterals  and  triangles.  The 
nodes  need  not  be  regularly  spaced;  the  configuration  of  the  elements  can  be  such  that  the  mesh 
closely  fits  an  irregularly  shaped  structure.  As  with  the  finite-difference  method,  there  is  a  trade¬ 
off  between  accuracy  (determined  by  the  number  of  elements,  i.e.  element  size)  and  cost  (deter¬ 
mined  by  the  amount  of  storage  and  the  number  of  calculations  required). 
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fillin'  14.  Forces  Q  ami  displacements  b  at  nodes  of  a  quadri¬ 
lateral  element  are  shown  for  the  eases  of  a)  SH  or  Love  waves 
and  h )  P  and  SV  t >r  Rayleigh  waves  ( from  l.ysmer  and  Drake  1972). 


An  assumption  inherent  to  the  finite-element  method  is  that  all  external  forces  and  forces  be¬ 
tween  elements  are  transmitted  through  nodal  points.  Each  nodal  point  has  one  or  two  force  com¬ 
ponents  acting  on  it,  depending  on  whether  the  wave  motion  considered  is  one-dimensional  (SH 
waves.  Love  waves)  or  two  dimensional  ( P  and  SF  waves,  Rayleigh  waves),  and  has  the  same  num¬ 
ber  of  displacement  components  (Fig.  14a,  b).  Three  requirements  that  must  be  satisfied  simul¬ 
taneously  are  those  of  equilibrium,  compatibility,  and  force-deflection  relationship,  i.e.  the  inter¬ 
nal  element  forces  acting  at  each  nodal  point  must  equilibrate  the  externally  applied  nodal  force, 
the  deformations  of  the  elements  must  be  such  that  the  elements  continue  to  meet  at  the  nodal 
points  in  the  loaded  condition,  and  the  internal  forces  and  internal  displacements  within  each 
element  must  be  related  as  required  by  its  individual  geometry  and  material  properties  (Clough 
1%5 1.  Displacements  of  points  within  an  element  are  defined  by  displacement  functions  in  terms 
of  the  displacements  at  nodes.  A  linear  element  is  one  for  which  the  internal  displacements  are 
linearly  interpolated  from  the  nodal  displacements.  Displacement  functions  are  chosen  so  that 
the  type  of  element  is  unchanged  as  a  disturbance  propagates  through  it  (i.e.  quadrilateral  elements 
remain  quadrilateral),  which  satisfies  the  displacement  compatibility  conditions  on  the  boundaries 
between  elements.  The  numerical  solution  is  obtained  either  for  the  case  of  a  forcing  function  or 
for  the  case  of  specified  initial  displacements  and  velocities. 

Displacements  of  nodal  points  are  related  to  forces  through  the  stiffness  matrix  [A'] .  If  the 
loading  is  static, 

[A|JM  =  )<?}•  (21) 

Mere  is  a  column  vector  which  represents  the  displacement  field;  the  size  of  |fi  |  is  determined 
by  the  number  of  nodal  points  (AT  and  by  the  degrees  of  freedom  at  each  nodal  point  ( }  5 }  is  of 
size  V  in  the  case  of  Love  or  SH  waves,  IN  in  the  case  of  P  and  SF  or  Rayleigh  waves).  The  column 
vector  jo |  represents  the  external  forces  applied  to  the  structure;  it  is  of  the  same  size  as  |5  |  . 
Once  [A  |  is  known,  eq  21  is  solved  for  tire  displacements  j  5  |  .  If  the  loading  is  dynamic  (time 
dependent ),  an  inertia  term  must  be  included  and 

[W|  l*"l  +  (AT  {*(  =  \V\  (22a) 

(Lysmcrand  Drake  1972).  If  a  damping  matrix  [C|  is  included. 


29 


(22b) 


[M]  |«"f  +  [C]  {S'}  +[A']{5}  =  \Q\ 

(Smith  1975).  Here  [A/]  is  a  mass  matrix.  The  prime  denotes  differentiation  with  respect  to  time. 
Once  [A/]  and  [A'  |  (and  [C] ,  if  included)  have  been  determined,  eq  22a,  b  may  be  solved  for  the 
displacements  j  6  j  .  When  the  excitation  of  the  structure  is  assumed  to  be  harmonic,  with  angular 
frequency  to,  steady -state  problems  can  be  solved  by  the  method  of  complex  response  (Lysmer 
and  Drake  1972).  Assume  that  j  P  |  is  a  time-independent,  complex  force  amplitude  and  juj  is 


a  time-independent,  complex  displacement  amplitude.  Let 

}(?}  =  {PjexpOW)  (23) 

and 

J  6  }  =  juj  exp(/co/).  (24) 

Substitution  of  eq  23  and  24  into  22b  results  in 

([A]  -to2  [A/])  ft/}  =  \P\  (25a) 

([A]  +  /«[C)  -  w2[A/])Jwj  =  jPj,  (25b) 

each  of  which  is  a  set  of  linear  equations  to  be  solved  for  displacement  amplitude  j  u  | .  This 


approach  is  correct  for  problems  involving  surface  waves  or  for  problems  involving  body  waves 
when  the  wavelength  is  greater  than  the  scale  size.  For  body  wave  problems  in  which  the  scale 
size  is  greater  than  the  wavelength,  eq  22a,  b  must  be  solved  in  the  time  domain  (Smith  1975). 
Lysmer  and  Drake  (1972)  describe  how  the  matrices  [M]  and  [A]  are  determined  for  the  solution 
in  the  frequency  domain.  The  procedure  involves  first  determining  the  mass  [A/*]  and  stiffness  [A’] 
matrices  for  each  element  and  then  assembling  the  element  matrices  to  form  [ M ]  and  [A] .  Because 
[Af*]  and  [A"]  are  determined  for  a  local  coordinate  system  convenient  for  the  configuration  of  the 
element,  they  must  be  transformed  into  a  global  coordinate  system  before  being  added  to  give  [M] 
and  [A] .  Smith  (1975)  uses  the  lumped  mass  matrix  which  is  preferred  for  problems  in  the  time  do¬ 
main. 

As  with  any  numerical  procedure,  the  accuracy  of  the  solution  obtained  with  the  finite-element 
method  must  be  evaluated.  Smith  (1975)  determined  that  for  the  case  of  a  plane  wave  propagat¬ 
ing  a  distance  of  1  km  through  linear,  rectangular  elements,  attenuation  is  negligible  and  disper¬ 
sion  is  =«  1%  for  a  mesh  size  of  10  elements  per  wavelength;  attenuation  and  dispersion  are  each 
3%  when  there  are  8  elements  per  wavelength.  These  results  are  dependent  on  the  time-stepping 
algorithm  and  on  the  length  of  the  time  step  Smith  used.  Eight  to  1 2  elements  and  6  to  1 0  ele¬ 
ments  per  wavelength  are  also  cited  in  the  literature.  Accuracy  is  also  affected  by  the  interpola¬ 
tion  used  to  define  displacements  within  an  element.  The  spacing  and  configuration  of  the  ele¬ 
ments  can  be  varied  in  accordance  wifh  the  resolution  needed  in  a  particular  portion  of  the  finite 
element  model.  Contamination  of  the  solution  by  waves  reflected  from  artificial  boundaries 
(sides  and  bottom  of  model)  is  avoided  through  the  use  of  nonreflecting  boundaries  or  by  termi¬ 
nating  the  solution  before  the  spurious  waves  reach  the  region  of  interest. 

Geller  et  al.  (1979)  determined  the  dynamic  finite-element  solutions  for  displacements  and 
velocities  observed  at  the  surface  due  to  movement  on  a  fault  modeled  as  a  dislocation.  Their  re¬ 
sults  are  compared  with  the  analytic  solution  in  Figures  15a  and  15b.  The  agreement  between  the 
analytic  solution  and  the  numerical  solution  is  worse  when  the  size  of  the  elements  is  doubled  be¬ 
cause  the  resolution  due  to  high  frequency  components  is  lost;  there  is  an  insufficient  number  of 
elements  per  wavelength  at  the  high  frequencies  when  the  element  size  is  increased.  The  agree¬ 
ment  between  the  analytic  and  the  numerical  solutions  is  better  for  displacements  than  for  veloc- 
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a.  Displacement  and  velocity  at  a  surface  location  0.5  km  above  and  5  km  off  the  fault. 
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b.  Displacement  and  velocity  at  the  surface  along  the  fault  strike.  The  fault  extends 
from  X  =  -5  km  to  X  =  5  km  at  a  depth  of  0.5  km. 


Figure  15.  Calculated  values  of  ground  motion  due  to  movement  along  a  right-lateral 
fault.  Results  are  given  for  the  analytic  solution  and  for  numerical  solutions  with  dif¬ 
ferent  element  sizes  (from  Geller  et  al.  1979). 
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ities  because  the  displacements  have  a  lower  predominant  frequency  and  are  thus  less  affected  by 
die  use  of  a  discrete  representation. 

Examples  of  problems  involving  wave  propagation  to  which  the  finite  element  is  well  suited  arc 
given  below. 

Lysmer  and  Drake  (1972)  apply  the  finite  element  method  to  the  study  of  the  propagation  of 
surface  waves  in  a  layered  medium  and  in  a  region  of  irregular  structure.  They  define  a  generalized 
Love  wave  to  be  any  free  motion  involving  harmonic  displacements  of  the  form  6  =  u(z)  exp  i(coi  - 
kx)  in  an  infinite  n-layered  structure  with  horizontal  interfaces.  Here  is  the  wave  number  and 
;/(r )  is  the  amplitude  function.  If  the  function  u(: )  is  assumed  to  vary  linearly  within  each  layer, 
then  a  discrete  representation  of  the  Love  wave  is 

{ S  }  =  a  |  u  }  exp  i(u>r  -  kx) , 

where  }  5 }  is  a  column  vector  containing  the  displacements  b-,j  =  1,2, ....  n  of  the  layer  interfaces, 
j  u  |  is  a  complex  amplitude  vector  with  the  components  u-rj  =1,2, ..., «,  and  a  is  a  mode  partici¬ 
pation  factor  which  indicates  the  contribution  of  each  mode  to  the  amplitude  of  the  displacement. 
The  assumption  of  linear  variation  of  the  amplitude  function  is  valid  if  the  thicknesses  of  the  model 
layers  are  small  compared  to  the  wavelengths  of  tire  shear  waves  in  the  layer  ;  this  requirement  may 
necessitate  introducing  more  layers  than  there  are  actual  structural  layers  present  in  the  model.  The 
equation  of  free  harmonic  motion  is 

(M]**  +  [G]  -cu2  [A/] )  {«}  =  jo}, 

where  the  matrices  [A )  and  [6'J  are  related  to  the  stiffness  of  the  layered  stmeture.  Boundary 
conditions  are  given  for  the  case  of  an  irregular  zone  between  two  semi-infinite  layered  zones.  The 
irregular  zone  has  surface  topography  and  layers  with  sloping  or  nonplanar  interfaces.  Vertical  in¬ 
terfaces  separate  the  irregular  zone  from  the  semi-infinite  layered  regions.  A  restriction  on  the 
thickness  of  layers  (relative  to  wavelength)  within  the  semi-infinite  zones  was  given  above;  an  addi¬ 
tional  restriction  is  that  the  interfaces  between  those  layers  must  coincide  with  nodal  points  on  the 
vertical  interfaces  with  the  irregular  zone.  Lysmer  and  Drake  (1972)  give  the  results  from  two 
problems,  propagation  of  Love  waves  across  an  idealized  alluvial  valley  and  propagation  of  Rayleigh 
waves  through  a  portion  of  tire  Central  Valley  of  California  and  the  Sierra  Nevada. 

Smith  ( 1 975 )  gives  the  results  of  a  finite-element  study  of  three  body  wave  propagation  prob¬ 
lems:  the  amplification  of  ground  motion  by  surface  topography  (mountain),  amplification  of 
ground  motion  due  to  a  semicylindrical  alluvial  valley,  and  ground  motion  due  to  a  deep  focus 
earthquake  (case  of  a  subducting  slab). 

Bolt  and  Smith  (1976)  use  the  finite-element  method  to  compute  synthetic  seismograms  for 
/’-waves  vertically  incident  on  a  massive  sulfide  ore  body  and  for  SH -waves  vertically  incident  on 
a  buried  salt  dome. 

McCowan  et  al.  (1977)  use  observed  ground  motions  from  the  San  Fernando,  California 
(9  Feb  1971 ),  earthquake  as  the  displacements  in  two-dimensional  finite-element  model  of  the 
fault.  They  first  consider  the  static  case  (eq  21)  and  solve  the  inverse  problem  for  the  force  vector. 
They  then  use  the  same  finite-element  model  for  the  dynamic  case  (eq  22b)  and  compute  synthet¬ 
ic  seismograms.  Finite-difference  equations  are  used  to  relate  acceleration  and  velocity  to  dis¬ 
placement  in  the  dynamic  problem. 

Melosh  and  Raefsky  (1981 )  utilize  the  technique  of  split  nodes  in  order  to  advantageously 
treat  tire  presence  of  faults  or  cracks  in  a  finite  element  model.  In  a  standard  finite-element 
scheme,  elements  that  share  a  node  experience  the  same  displacement  at  that  node;  the  distinguish¬ 
ing  characteristic  of  a  split  node  is  that  its  displacement  varies  with  the  element  under  considera¬ 
tion.  In  this  way,  a  split  node  located  on  a  fault  has  two  displacement  values,  the  first  (n+)  ex- 
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perienced  by  the  element  on  one  side  of  the  fault,  the  second  (;/')  experienced  by  the  element  on 
the  opposite  side  of  the  fault.  Melosh  and  Ruefsky  state  that  the  use  of  split  nodes  does  not  com¬ 
promise  the  accuracy  or  stability  of  the  computation. 


SECTION  8.  HYBRID  METHODS 

The  characteristics  of  a  problem  may  determine  that  it  is  advantageous  to  employ  a  method  of 
generating  synthetic  seismograms  that  is  based  upon  ray  theory  rather  titan  mode  theory,  or  vice 
versa.  Table  2  (from  Urick  1983,  p.  122)  compares  normal  mode  and  ray  theories.  A  normal 
mode  approach  handles  initial  conditions  well  and  is  useful  at  low  frequencies  (below  500-1000 
Hz  depending  on  the  model  used),  while  a  method  based  upon  ray  theory  handles  boundary  condi¬ 
tions  well  and  is  appropriate  for  programs  involving  high  frequencies  and  long  propagation  dis¬ 
tances;  although  mode  theory  is  valid  at  all  frequencies,  its  use  at  high  frequencies  effectively  is 
restricted  to  short  propagation  distances.  Hybrid  methods  that  have  features  of  both  ray  theory 
and  normal  mode  theory  may  have  lower  computational  cost  yet  retain  sufficient  accuracy. 
RAYMODE,  a  computer  program  written  for  underwater  acoustics,  is  a  hybrid  method  applicable 
to  studies  of  compressional  wave  propagation  at  a  range  of  frequencies  (low  to  high)  in  a  medium 
characterized  by  depth-dependent  velocity  and  a  flat  lower  boundary  at  which  there  may  be  loss 
of  energy. 

Another  hybrid  ray-mode  formulation  developed  for  studies  of  wave  propagation  in  wave¬ 
guides  and  ducts  ( Felsen  and  Ishihara  1979,  Felsen  1 98 1 )  has  been  used  to  construct  synthetic 
seismograms  for  the  case  of  SH- waves  in  a  layer  over  a  half-space  (Kamel  and  Felsen  1981).  This 
method  employs  a  truncated  ray  series,  each  ray  characterized  by  its  takeoff  angle  at  the  source, 
and  a  truncated  mode  series,  each  mode  characterized  by  the  angle  of  incidence  of  its  equivalent 
plane  waves  at  the  layer  boundaries;  together,  the  two  series  constitute  a  complete  spectrum  of 
angles  while  avoiding  problems  such  as  caustics  by  the  appropriate  selection  of  rays  to  be  incorpor¬ 
ated  in  the  solution.  Because  the  total  number  of  generalized  rays  and  normal  modes  contained  in 
this  hybrid  formulation  generally  is  less  than  the  number  of  rays  or  of  modes  required  when  either 
method  is  employed  independently,  the  computational  cost  of  the  hybrid  method  is  lower.  This 
hybrid  method,  developed  for  the  case  of  high  frequencies,  includes  asymptotic  approximations 
and  results  in  a  solution  of  the  wave  equation  comprised  of  generalized  rays,  normal  modes,  and  a 
remainder;  Felsen  (1981 )  notes  that  if  the  asymptotic  approximations  are  not  used,  the  hybrid 
formulation  is  exact.  Synthetic  seismograms  generated  by  four  methods  were  compared.  The  hy¬ 
brid  formulation  consisted  of  two  generalized  rays  (direct  and  singly  reflected),  two  trapped  modes, 


Table  2.  Comparison  of  normal-mode  and  ray  theories  (from  Urick  1983). 


Normal'nwde  theory  Ray  theory 


Gives  a  formally  complete  solution. 

Solution  is  difficult  to  interpret. 

Cannot  easily  handle  real  boundary 
conditions. 

Source  function  easily  inserted. 

Requires  a  computer  program,  except 
in  limiting  cases  when  analytic  an¬ 
swers  exist,  and  presents  computa¬ 
tional  difficulties  in  all  but  simplest 
boundary  conditions. 

Valid  at  all  frequencies,  but  practically 
is  useful  for  low  frequencies  (few 
modes). 


Does  not  handle  diffraction  problems, 
e.g.  sound  in  a  shadow. 

Rays  are  easily  drawn.  Sound  distribu¬ 
tion  is  easily  visualized. 

Real  boundary  conditions  are  inserted 
easily,  e.g.  a  sloping  bottom. 

Is  independent  of  the  source. 

Rays  can  be  drawn  by  hand  using  Snell’s 
law.  However,  a  ray-trace  computer 
program  is  normally  used. 


Valid  only  at  high  frequencies  if  (1) 
radius  of  ray  curvature  is  >  \  or  (2) 
sound  velocity  does  not  change  much 
in  a  i.e.  (AV{V)/Az  «  |/\. _ 


and  a  remainder.  The  method  of  generalized  rays  (Cagniard-de  Hoop)  incorporated  14  rays.  Seismo¬ 
grams  generated  by  these  two  methods  were  indistinguishable.  The  seismogram  obtained  by  a  solu¬ 
tion  equivalent  to  the  hybrid  formulation  without  the  remainder  term  agreed  with  that  of  the  hy¬ 
brid  method  for  early  observation  times  but  became  in  error  at  a  time  corresponding  to  the  arrival 
of  the  singly  reflected  cay;  from  then  on,  the  discrepancy  between  t\e  two  seismograms  decreased 
with  time.  The  fourth  seismogram,  based  upon  the  two  trapped  modes  incorporated  in  the  hybrid 
solution,  was  initially  in  error  but  approached  agreement  with  that  of  the  hybrid  method  with 
time. 

A  second  type  of  hyorid  method  is  one  which  combines  analytical  and  numerical  methods  of 
solving  the  wave  equation.  Shtivelman  (1984)  uses  the  finite  difference  method  and  a  form  of 
asymptotic  ray  theory  to  compute  the  wave  field  in  media  having  irregular  structures.  The  finite 
difference  method  is  used  in  the  vicinity  of  a  scatterer  or  complex  structure  because  the  numerical 
method  generally  handles  wave  propagation  in  heterogeneous  media  well,  whereas  the  ray  method 
is  unsuitable  when  scattering  from  heterogeneities  or  wave  propagation  through  thin  beds  is  in¬ 
volved.  The  ray  method  is  used  to  compute  the  wave  field  in  homogeneous  regions  or  in  regions 
where  layers  are  thick  relative  to  wavelength  because  in  such  cases  the  computational  cost  and 
storage  requirements  of  the  finite  difference  method,  together  with  the  possibility  of  numeric  dis¬ 
persion  causing  loss  of  accuracy,  cause  it  to  be  the  less  suitable  method.  The  transition  between 
methods  occurs  at  an  arbitrary  contour  along  which  the  wave  field  of  the  first  method  employed 
is  used  to  define  source  terms  for  the  second  method.  Shtivelman  gives  examples  of  the  applica¬ 
tion  of  this  hybrid  method  to  the  case  of  SH  plane  waves  impinging  on  an  interface  that  is  displaced 
by  a  right  angle  step,  an  inclined  step,  a  right  angle  projection,  a  right  angle  recess,  or  a  thrust  fault. 


SECTION  9.  CONCLUSION 


Several  methods  of  generating  synthetic  seismograms  have  been  reviewed.  The  selection  of  a 
method  for  use  is  based  upon  considerations  of  the  physical  aspects  of  a  problem -complexity  of 
the  geologic  model,  wave  types  to  be  included  in  the  seismogram  and  of  computational  cost.  The 
type  of  problem  for  which  various  methods  are  appropriate  is  indicated  by  Table  3,  in  which  speci¬ 
fic  computer  codes  are  evaluated  by  Braile  (1982),  and  by  Table  4. 

Table  3.  Summary  of  seismic  interpretation  techniques  (from  Braile  1982). 


Method 

Characteristics 

Limitations 

Travel-time  Analysis 

Geometrical  Ray  theory 
calculation  of  travel-times 
of  reflected  and  refracted 
phases  (Program  TXCURV). 

Travel-times  for  head 
waves  and  reflections 
for  plane  or  dipping 
layered  models;  exact. 

Model  geometry; 
travel-time  only. 

Ray  Tracing  (Program 

RAY2D). 

Two-dimensional  models. 

Travel-times  only; 
limited  number  of 
phases;  ignores  PS 
conversions,  guided 
waves  and  surface 
waves. 

Synthetic  Seismogram  Modeling 

Reflectivity  (Program 

SYNCAL). 

Wave-theoretical  solution 
to  response  of  a  layered 

Restricted  to  1-D 
models. 

half-space.  Includes  P,  SK, 
guided  phases  and  surface 
waves.  Exact.  Includes  Q. 
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Table  3  (cont’d). 


Table  4  (cont’d).  Characteristics  of  methods  of  generating  synthetic  seismograms. 

Method _  Characteristics 

Surface  wave  seismology  frequencies  that  are  too  low  to  meet  high  frequency  approxi- 

(cont’d)  mation  of  ray/wave  theory  yet  are  high  enough  that  required 

number  of  normal  modes  is  prohibitive. 

Normal  modes  Analytic.  Plane-layered  waveguide  handled  readily.  Valid 

at  all  frequencies  but  practical  at  low  frequencies  (smaller 
number  of  modes  and  so  lower  cost).  Computational  time 
reduced  by  employing  Fast  Field  Program  which  can  in¬ 
corporate  both  normal  and  leaky  modes.  Costly,  perhaps 
inefficient,  in  case  of  complex  structure  because  normal 
modes  must  be  recomputed  whenever  inhomogeneity  in 
model  is  encountered. 

Numerical.  Displacement  at  grid  points.  Useful  in  near- 
field  or  with  complex  structure.  Body  waves  or  surface 
waves.  Size  of  grid  spacing  affects  stability  and  accuracy 
of  results.  Truncation  error  due  to  approximation  to 
wave  equation.  Aliasing,  numeric  dispersion,  and  spurious 
reflected  waves  from  artificial  boundaries  are  potential 
problems.  Costly  due  to  number  of  computations  and 
storage  requirements. 

Finite  elements  Numerical.  Displacements  between  grid  points  known  by 

interpolation.  Useful  with  irregularly  shaped  structure. 

Body  waves  or  surface  waves.  Accuracy  dependent  on 
lumber  of  elements  and  method  of  interpolation.  Alias¬ 
ing,  numeric  dispersion,  and  spurious  reflected  waves 
from  artificial  boundaries  are  potential  problems.  Costly 
_ due  to  number  of  computations  and  storage  requirements. 


By  presenting  the  characteristic  features  of  each  method,  this  report  may  overemphasize  differ¬ 
ences  and  unintentionally  obscure  similarities  among  the  methods.  It  should  be  realized  that  the 
accuracy  of  synthetic  seismograms  generated  by  any  of  these  methods  is  highly  frequency  depend¬ 
ent.  Frequency  (f)  and  equivalently  wavelength  (X  =  V/f)  serve  as  scale  factors  for  a  problem.  The 
validity  of  the  synthetic  seismogram  obtained  is  constrained  by  grid  spacing  relative  to  wavelength 
with  numerical  methods  and  by  layer  thickness  relative  to  wavelength  or  change  in  velocity  gradi¬ 
ent  over  a  wavelength  in  analytic  approaches  to  solving  the  wave  equation.  The  condition  for  the 
valid  use  of  asymptotic  approximations  to  replace  exact  quantities  in  the  analytic  solutions  is  gen¬ 
erally  kr  -  2rrr/X  »  1 ,  where  r  is  the  range  from  the  source.  In  the  far  field,  r  is  sufficiently  large 
(at  least  10  X)  that  the  condition  is  satisfied.  The  actual  propagation  distance  required  for  far -field 
conditions  to  prevail  is  inversely  proportional  to  frequency. 

An  understanding  of  the  limitations  and  approximations  inherent  to  the  methods  of  generating 
synthetic  seismograms,  together  with  the  selection  of  a  method  appropriate  to  a  particular  problem, 
is  necessary  if  synthetic  seismograms  are  to  be  used  wisely  and  advantageously. 
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